Skip to content

Fix legacy vertexing: drop ill-conditioned tracks and skip unneeded momenta-at-vertex computation #525

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
kjvbrt wants to merge 5 commits into
base: pre-edm4hep1
Choose a base branch
from
Open

Fix legacy vertexing: drop ill-conditioned tracks and skip unneeded momenta-at-vertex computation #525

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c3f1e86
Drop tracks with non positive-definite covariance before vertex fitting
kjvbrt Jul 20, 2026
ddda32c
Applyied clang formatting
kjvbrt Jul 20, 2026
ac64cb5
Make momenta-at-vertex computation in VertexFitter_Tk optional
kjvbrt Jul 21, 2026
8209073
Thread ComputeMomentaAtVertex through VertexFitter() and myUtils::get…
kjvbrt Jul 21, 2026
9c83e46
Run clang-format
kjvbrt Jul 21, 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
53 changes: 33 additions & 20 deletions analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,26 +32,39 @@ namespace FCCAnalyses{

namespace VertexFitterSimple{

/// Vertex (code from Franco Bedeschi): passing the recoparticles. Units for the beamspot constraint: mum
VertexingUtils::FCCAnalysesVertex VertexFitter( int Primary,
ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> recoparticles,
ROOT::VecOps::RVec<edm4hep::TrackState> alltracks,
bool BeamSpotConstraint = false,
double sigmax=0., double sigmay=0., double sigmaz=0.,
double bsc_x=0., double bsc_y=0., double bsc_z=0. ) ;


/// Vertex (code from Franco Bedeschi): passing the tracks. Units for the beamspot constraint: mum
VertexingUtils::FCCAnalysesVertex VertexFitter_Tk( int Primary, ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
bool BeamSpotConstraint = false,
double sigmax=0., double sigmay=0., double sigmaz=0.,
double bsc_x=0., double bsc_y=0., double bsc_z=0. ) ;

VertexingUtils::FCCAnalysesVertex VertexFitter_Tk( int Primary, ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
const ROOT::VecOps::RVec<edm4hep::TrackState>& alltracks,
bool BeamSpotConstraint = false,
double sigmax=0., double sigmay=0., double sigmaz=0.,
double bsc_x=0., double bsc_y=0., double bsc_z=0. ) ;
/// Vertex (code from Franco Bedeschi): passing the recoparticles. Units for the
/// beamspot constraint: mum. See VertexFitter_Tk() below for what
/// ComputeMomentaAtVertex controls.
VertexingUtils::FCCAnalysesVertex VertexFitter(
int Primary,
ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> recoparticles,
ROOT::VecOps::RVec<edm4hep::TrackState> alltracks,
bool BeamSpotConstraint = false, double sigmax = 0., double sigmay = 0.,
double sigmaz = 0., double bsc_x = 0., double bsc_y = 0., double bsc_z = 0.,
bool ComputeMomentaAtVertex = true);

/// Vertex (code from Franco Bedeschi): passing the tracks. Units for the
/// beamspot constraint: mum. ComputeMomentaAtVertex controls whether
/// updated_track_momentum_at_vertex is filled: it requires building a Delphes
/// VertexMore object, whose cost scales as O(ntracks^3) and dominates the total
/// runtime for events with many tracks. Set to false when only the vertex
/// position/chi2 is needed (e.g. FCCAnalysesVertex::vertex via
/// VertexingUtils::get_VertexData) -- see
/// https://github.com/HEP-FCC/FCCAnalyses/issues/378.
VertexingUtils::FCCAnalysesVertex
VertexFitter_Tk(int Primary, ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
bool BeamSpotConstraint = false, double sigmax = 0.,
double sigmay = 0., double sigmaz = 0., double bsc_x = 0.,
double bsc_y = 0., double bsc_z = 0.,
bool ComputeMomentaAtVertex = true);

VertexingUtils::FCCAnalysesVertex
VertexFitter_Tk(int Primary, ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
const ROOT::VecOps::RVec<edm4hep::TrackState> &alltracks,
bool BeamSpotConstraint = false, double sigmax = 0.,
double sigmay = 0., double sigmaz = 0., double bsc_x = 0.,
double bsc_y = 0., double bsc_z = 0.,
bool ComputeMomentaAtVertex = true);

/// Return the tracks that are flagged as coming from the primary vertex
ROOT::VecOps::RVec<edm4hep::TrackState> get_PrimaryTracks( ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
Expand Down
17 changes: 12 additions & 5 deletions analyzers/dataframe/FCCAnalyses/myUtils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,11 +107,18 @@ namespace myUtils{
ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertex>
merge_VertexObjet(ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertex> in);

ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertex> get_VertexObject(ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertexMC> mcver,
ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> reco,
ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
ROOT::VecOps::RVec<int> recin,
ROOT::VecOps::RVec<int> mcin);
/// ComputeMomentaAtVertex controls whether updated_track_momentum_at_vertex
/// is filled for each returned vertex: it requires building a Delphes
/// VertexMore object, whose cost scales as O(ntracks^3) and dominates the
/// total runtime for events with many tracks. Set to false when only the
/// vertex position/chi2/primary flag is needed (e.g. hasPV()) -- see
/// https://github.com/HEP-FCC/FCCAnalyses/issues/378.
ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertex> get_VertexObject(
ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertexMC> mcver,
ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> reco,
ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
ROOT::VecOps::RVec<int> recin, ROOT::VecOps::RVec<int> mcin,
bool ComputeMomentaAtVertex = true);

ROOT::VecOps::RVec<float> get_Vertex_mass(ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertex> vertex,
ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> reco);
Expand Down
88 changes: 70 additions & 18 deletions analyzers/dataframe/src/VertexFitterSimple.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@

#include <fcntl.h>

//#include "TrkUtil.h" // from delphes
#include "TrackCovariance/TrkUtil.h" // from Delphes

namespace FCCAnalyses {

Expand All @@ -35,12 +35,44 @@ void resume_stdout(int fd) {

// -----------------------------------------------------------------------------

// Drop tracks whose covariance matrix is not positive-definite: VertexFit
// (Delphes) silently produces a NaN vertex when such a track is included,
// instead of failing loudly, on top of flooding stdout while it tries (and
// fails) to invert the matrix. See
// https://github.com/HEP-FCC/FCCAnalyses/issues/378
//
// Printed via std::cerr rather than std::cout since callers wrap the fit
// itself in supress_stdout()/resume_stdout() to silence Delphes' own
// printf() spam; stderr (fd 2) is untouched by that redirection, so this
// stays visible without reintroducing the spam.
static ROOT::VecOps::RVec<edm4hep::TrackState>
drop_track_with_bad_covariance(ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
bool Units_mm) {
ROOT::VecOps::RVec<edm4hep::TrackState> good_tracks;
for (auto &t : tracks) {
TMatrixDSym Cov = VertexingUtils::get_trackCov(t, Units_mm);
if (TrkUtil::CheckPosDef(Cov)) {
good_tracks.push_back(t);
continue;
}
std::cerr
<< "VertexFitterSimple: dropping track with non positive-definite "
"covariance matrix -- track parameters: D0="
<< t.D0 << " phi=" << t.phi << " omega=" << t.omega << " Z0=" << t.Z0
<< " tanLambda=" << t.tanLambda
<< " -- covariance diagonal (D0,phi,omega,Z0,tanLambda)=(" << Cov(0, 0)
<< ", " << Cov(1, 1) << ", " << Cov(2, 2) << ", " << Cov(3, 3) << ", "
<< Cov(4, 4) << ")" << std::endl;
}
return good_tracks;
}

VertexingUtils::FCCAnalysesVertex VertexFitter(
int Primary,
ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> recoparticles,
ROOT::VecOps::RVec<edm4hep::TrackState> thetracks, bool BeamSpotConstraint,
double bsc_sigmax, double bsc_sigmay, double bsc_sigmaz, double bsc_x,
double bsc_y, double bsc_z) {
double bsc_y, double bsc_z, bool ComputeMomentaAtVertex) {

// input = a collection of recoparticles (in case one want to make
// associations to RecoParticles ?) and thetracks = the collection of all
Expand All @@ -56,9 +88,9 @@ VertexingUtils::FCCAnalysesVertex VertexFitter(

// FCCAnalysesVertex thevertex = VertexFitter_Tk( Primary, tracks, thetracks)
// ;
thevertex =
VertexFitter_Tk(Primary, tracks, thetracks, BeamSpotConstraint,
bsc_sigmax, bsc_sigmay, bsc_sigmaz, bsc_x, bsc_y, bsc_z);
thevertex = VertexFitter_Tk(Primary, tracks, thetracks, BeamSpotConstraint,
bsc_sigmax, bsc_sigmay, bsc_sigmaz, bsc_x, bsc_y,
bsc_z, ComputeMomentaAtVertex);

// fill the indices of the tracks
ROOT::VecOps::RVec<int> reco_ind;
Expand All @@ -84,11 +116,13 @@ VertexingUtils::FCCAnalysesVertex VertexFitter(
VertexingUtils::FCCAnalysesVertex
VertexFitter_Tk(int Primary, ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
bool BeamSpotConstraint, double bsc_sigmax, double bsc_sigmay,
double bsc_sigmaz, double bsc_x, double bsc_y, double bsc_z) {
double bsc_sigmaz, double bsc_x, double bsc_y, double bsc_z,
bool ComputeMomentaAtVertex) {

ROOT::VecOps::RVec<edm4hep::TrackState> dummy;
return VertexFitter_Tk(Primary, tracks, dummy, BeamSpotConstraint, bsc_sigmax,
bsc_sigmay, bsc_sigmaz, bsc_x, bsc_y, bsc_z);
bsc_sigmay, bsc_sigmaz, bsc_x, bsc_y, bsc_z,
ComputeMomentaAtVertex);
}

// ---------------------------------------------------------------------------------------------------------------------------
Expand All @@ -97,13 +131,18 @@ VertexingUtils::FCCAnalysesVertex
VertexFitter_Tk(int Primary, ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
const ROOT::VecOps::RVec<edm4hep::TrackState> &alltracks,
bool BeamSpotConstraint, double bsc_sigmax, double bsc_sigmay,
double bsc_sigmaz, double bsc_x, double bsc_y, double bsc_z) {
double bsc_sigmaz, double bsc_x, double bsc_y, double bsc_z,
bool ComputeMomentaAtVertex) {
// Suppressing printf() output from TMatrixBase:
// https://github.com/root-project/root/blob/722eb4652bfc79149df00c8b0e92d0837caf054c/math/matrix/src/TMatrixTBase.cxx#L662
// The solution found here:
// https://stackoverflow.com/questions/46728680/how-to-temporarily-suppress-output-from-printf
int fd = supress_stdout();

bool Units_mm = true;

tracks = drop_track_with_bad_covariance(tracks, Units_mm);

// Units for the beam-spot : mum
// See
// https://github.com/HEP-FCC/FCCeePhysicsPerformance/tree/master/General#generating-events-under-realistic-fcc-ee-environment-conditions
Expand Down Expand Up @@ -151,8 +190,6 @@ VertexFitter_Tk(int Primary, ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
TVectorD **trkPar = new TVectorD *[Ntr];
TMatrixDSym **trkCov = new TMatrixDSym *[Ntr];

bool Units_mm = true;

for (Int_t i = 0; i < Ntr; i++) {
edm4hep::TrackState t = tracks[i];
TVectorD par = VertexingUtils::get_trackParam(t, Units_mm);
Expand Down Expand Up @@ -213,20 +250,25 @@ VertexFitter_Tk(int Primary, ROOT::VecOps::RVec<edm4hep::TrackState> tracks,

TheVertex.vertex = result;

// Use VertexMore to retrieve more information :
VertexMore theVertexMore(&theVertexFit, Units_mm);

for (Int_t i = 0; i < Ntr; i++) {

TVectorD updated_par =
theVertexFit.GetNewPar(i); // updated track parameters
TVectorD updated_par_edm4hep =
VertexingUtils::Delphes2Edm4hep_TrackParam(updated_par, Units_mm);
updated_track_parameters.push_back(updated_par_edm4hep);
}

// Momenta of the tracks at the vertex:
TVector3 ptrack_at_vertex = theVertexMore.GetMomentum(i);
updated_track_momentum_at_vertex.push_back(ptrack_at_vertex);
// VertexMore computes the joint covariance of the vertex position and
// every track's momentum at the vertex, which costs O(Ntr^3) -- this
// dominates the total runtime for events with many tracks (see
// https://github.com/HEP-FCC/FCCAnalyses/issues/378). Skip it whenever
// the caller doesn't need updated_track_momentum_at_vertex.
if (ComputeMomentaAtVertex) {
VertexMore theVertexMore(&theVertexFit, Units_mm);
for (Int_t i = 0; i < Ntr; i++) {
TVector3 ptrack_at_vertex = theVertexMore.GetMomentum(i);
updated_track_momentum_at_vertex.push_back(ptrack_at_vertex);
}
}

TheVertex.updated_track_parameters = updated_track_parameters;
Expand Down Expand Up @@ -274,10 +316,18 @@ get_PrimaryTracks(ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
<< tracks.size() << std::endl;
}

// Suppressing printf() output from TMatrixBase / Delphes' VertexFit, see
// VertexFitter_Tk() above for details.
int fd = supress_stdout();

tracks = drop_track_with_bad_covariance(tracks, false);

ROOT::VecOps::RVec<edm4hep::TrackState> seltracks = tracks;

if (seltracks.size() <= 1)
if (seltracks.size() <= 1) {
resume_stdout(fd);
return seltracks;
}

int Ntr = tracks.size();

Expand Down Expand Up @@ -344,6 +394,8 @@ get_PrimaryTracks(ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
delete[] trkPar;
delete[] trkCov;

resume_stdout(fd);

return seltracks;
}

Expand Down
19 changes: 12 additions & 7 deletions analyzers/dataframe/src/myUtils.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -239,13 +239,12 @@ ROOT::VecOps::RVec<int> get_Vertex_indMC(ROOT::VecOps::RVec<VertexingUtils::FCCA

}


ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertex>
get_VertexObject(ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertexMC> mcver,
ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> reco,
ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
ROOT::VecOps::RVec<int> recin,
ROOT::VecOps::RVec<int> mcin){
ROOT::VecOps::RVec<edm4hep::ReconstructedParticleData> reco,
ROOT::VecOps::RVec<edm4hep::TrackState> tracks,
ROOT::VecOps::RVec<int> recin, ROOT::VecOps::RVec<int> mcin,
bool ComputeMomentaAtVertex) {

ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertex> result;

Expand All @@ -272,9 +271,15 @@ get_VertexObject(ROOT::VecOps::RVec<VertexingUtils::FCCAnalysesVertexMC> mcver,
if (recoparticles.size()<2)continue;

VertexingUtils::FCCAnalysesVertex TheVertex;
if (v==0) TheVertex = VertexFitterSimple::VertexFitter(1,recoparticles, tracks, true, 4.5, 20e-3, 300 );
if (v == 0)
TheVertex = VertexFitterSimple::VertexFitter(
1, recoparticles, tracks, true, 4.5, 20e-3, 300, 0., 0., 0.,
ComputeMomentaAtVertex);
//if (v==0) TheVertex = VertexFitterSimple::VertexFitter(1,recoparticles, tracks);
else TheVertex = VertexFitterSimple::VertexFitter(0,recoparticles, tracks);
else
TheVertex = VertexFitterSimple::VertexFitter(0, recoparticles, tracks,
false, 0., 0., 0., 0., 0.,
0., ComputeMomentaAtVertex);
if (std::isnan(TheVertex.vertex.chi2))continue;
TheVertex.mc_ind=v;
result.push_back(TheVertex);
Expand Down
14 changes: 10 additions & 4 deletions examples/FCCee/tutorials/vertexing/analysis_primary_vertex.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,13 +27,18 @@ def analysers(df):
# std::array<float, 6> covMatrix{}; ///< covariance matrix of the position (stored as lower triangle matrix, i.e. cov(xx),cov(y,x),cov(z,x),cov(y,y),... )
# - ROOT::VecOps::RVec<float> reco_chi2 : the contribution to the chi2 of all tracks used in the fit
# - ROOT::VecOps::RVec< TVector3 > updated_track_momentum_at_vertex : the post-fit (px, py, pz ) of the tracks, at the vertex (and not at their d.c.a.)
.Define("VertexObject_allTracks", "VertexFitterSimple::VertexFitter_Tk ( 1, EFlowTrack_1, true, 4.5, 20e-3, 300)")
# (left empty below via ComputeMomentaAtVertex=false, since only the vertex position/chi2 is used
# in this tutorial -- computing it requires a Delphes VertexMore object whose cost scales as
# O(ntracks^3) and dominates the runtime for events with many tracks, see
# https://github.com/HEP-FCC/FCCAnalyses/issues/378)
.Define("VertexObject_allTracks", "VertexFitterSimple::VertexFitter_Tk ( 1, EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0., false)")

# EFlowTrack_1 is the collection of all tracks (the fitting method can of course be applied to a subset of tracks (see later)).
# "true" means that a beam-spot constraint is applied. Default is no BSC. Following args are the BS size and position, in mum :
# bool BeamSpotConstraint = false,
# double sigmax=0., double sigmay=0., double sigmaz=0.,
# double bsc_x=0., double bsc_y=0., double bsc_z=0. ) ;
# double bsc_x=0., double bsc_y=0., double bsc_z=0.,
# bool ComputeMomentaAtVertex = true ) ;

# This returns the edm4hep::VertexData :
.Define("Vertex_allTracks", "VertexingUtils::get_VertexData( VertexObject_allTracks )") # primary vertex, in mm
Expand All @@ -42,8 +47,9 @@ def analysers(df):
# Below, we determine the "primary tracks" using an iterative algorithm - cf LCFI+.
.Define("RecoedPrimaryTracks", "VertexFitterSimple::get_PrimaryTracks( EFlowTrack_1, true, 4.5, 20e-3, 300, 0., 0., 0.)")

# Now we run again the vertex fit, but only on the primary tracks :
.Define("PrimaryVertexObject", "VertexFitterSimple::VertexFitter_Tk ( 1, RecoedPrimaryTracks, true, 4.5, 20e-3, 300) ")
# Now we run again the vertex fit, but only on the primary tracks
# (again skipping the momenta-at-vertex computation, unused here):
.Define("PrimaryVertexObject", "VertexFitterSimple::VertexFitter_Tk ( 1, RecoedPrimaryTracks, true, 4.5, 20e-3, 300, 0., 0., 0., false) ")
.Define("PrimaryVertex", "VertexingUtils::get_VertexData( PrimaryVertexObject )")

# It is often useful to retrieve the secondary (i.e. non-primary) tracks, for example to search for secondary vertices.
Expand Down
Loading