diff --git a/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h b/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h index 257a5f9c488..a59ce5f4237 100644 --- a/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h +++ b/analyzers/dataframe/FCCAnalyses/VertexFitterSimple.h @@ -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 recoparticles, - ROOT::VecOps::RVec 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 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 tracks, - const ROOT::VecOps::RVec& 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 recoparticles, + ROOT::VecOps::RVec 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 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 tracks, + const ROOT::VecOps::RVec &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 get_PrimaryTracks( ROOT::VecOps::RVec tracks, diff --git a/analyzers/dataframe/FCCAnalyses/myUtils.h b/analyzers/dataframe/FCCAnalyses/myUtils.h index 2b8f3a6ded2..f60f1b2c7b9 100644 --- a/analyzers/dataframe/FCCAnalyses/myUtils.h +++ b/analyzers/dataframe/FCCAnalyses/myUtils.h @@ -107,11 +107,18 @@ namespace myUtils{ ROOT::VecOps::RVec merge_VertexObjet(ROOT::VecOps::RVec in); - ROOT::VecOps::RVec get_VertexObject(ROOT::VecOps::RVec mcver, - ROOT::VecOps::RVec reco, - ROOT::VecOps::RVec tracks, - ROOT::VecOps::RVec recin, - ROOT::VecOps::RVec 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 get_VertexObject( + ROOT::VecOps::RVec mcver, + ROOT::VecOps::RVec reco, + ROOT::VecOps::RVec tracks, + ROOT::VecOps::RVec recin, ROOT::VecOps::RVec mcin, + bool ComputeMomentaAtVertex = true); ROOT::VecOps::RVec get_Vertex_mass(ROOT::VecOps::RVec vertex, ROOT::VecOps::RVec reco); diff --git a/analyzers/dataframe/src/VertexFitterSimple.cc b/analyzers/dataframe/src/VertexFitterSimple.cc index 3249ae62c4a..8e728ddfcb8 100644 --- a/analyzers/dataframe/src/VertexFitterSimple.cc +++ b/analyzers/dataframe/src/VertexFitterSimple.cc @@ -8,7 +8,7 @@ #include -//#include "TrkUtil.h" // from delphes +#include "TrackCovariance/TrkUtil.h" // from Delphes namespace FCCAnalyses { @@ -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 +drop_track_with_bad_covariance(ROOT::VecOps::RVec tracks, + bool Units_mm) { + ROOT::VecOps::RVec 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 recoparticles, ROOT::VecOps::RVec 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 @@ -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 reco_ind; @@ -84,11 +116,13 @@ VertexingUtils::FCCAnalysesVertex VertexFitter( VertexingUtils::FCCAnalysesVertex VertexFitter_Tk(int Primary, ROOT::VecOps::RVec 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 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); } // --------------------------------------------------------------------------------------------------------------------------- @@ -97,13 +131,18 @@ VertexingUtils::FCCAnalysesVertex VertexFitter_Tk(int Primary, ROOT::VecOps::RVec tracks, const ROOT::VecOps::RVec &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 @@ -151,8 +190,6 @@ VertexFitter_Tk(int Primary, ROOT::VecOps::RVec 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); @@ -213,20 +250,25 @@ VertexFitter_Tk(int Primary, ROOT::VecOps::RVec 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; @@ -274,10 +316,18 @@ get_PrimaryTracks(ROOT::VecOps::RVec 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 seltracks = tracks; - if (seltracks.size() <= 1) + if (seltracks.size() <= 1) { + resume_stdout(fd); return seltracks; + } int Ntr = tracks.size(); @@ -344,6 +394,8 @@ get_PrimaryTracks(ROOT::VecOps::RVec tracks, delete[] trkPar; delete[] trkCov; + resume_stdout(fd); + return seltracks; } diff --git a/analyzers/dataframe/src/myUtils.cc b/analyzers/dataframe/src/myUtils.cc index e348dd41e83..f6776057ce6 100644 --- a/analyzers/dataframe/src/myUtils.cc +++ b/analyzers/dataframe/src/myUtils.cc @@ -239,13 +239,12 @@ ROOT::VecOps::RVec get_Vertex_indMC(ROOT::VecOps::RVec get_VertexObject(ROOT::VecOps::RVec mcver, - ROOT::VecOps::RVec reco, - ROOT::VecOps::RVec tracks, - ROOT::VecOps::RVec recin, - ROOT::VecOps::RVec mcin){ + ROOT::VecOps::RVec reco, + ROOT::VecOps::RVec tracks, + ROOT::VecOps::RVec recin, ROOT::VecOps::RVec mcin, + bool ComputeMomentaAtVertex) { ROOT::VecOps::RVec result; @@ -272,9 +271,15 @@ get_VertexObject(ROOT::VecOps::RVec 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); diff --git a/examples/FCCee/tutorials/vertexing/analysis_primary_vertex.py b/examples/FCCee/tutorials/vertexing/analysis_primary_vertex.py index e95fcd58d39..623e3a54ba6 100644 --- a/examples/FCCee/tutorials/vertexing/analysis_primary_vertex.py +++ b/examples/FCCee/tutorials/vertexing/analysis_primary_vertex.py @@ -27,13 +27,18 @@ def analysers(df): # std::array 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 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 @@ -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.