Fix minimize function in Measurement.cpp

- remove dead code
- fix and document minimize function in measurement.cpp
- make documentation a bit clearer. Thanks to @ldufour for feedback!

Tr/TrackFitEvent/src/Measurement.cpp

@@ -75,21 +75,58 @@ namespace {
 
 LHCb::Minimize1DResult LHCb::minimize( const LHCb::Measurement& m, const LHCb::StateZTraj<double>& refTraj,
                                        double zState ) {
-  // assume minimum is close to zState
-  auto       dist   = refTraj.position( zState ) - m.trajectory().position( 0 );
-  const auto d0     = m.trajectory().direction( 0 );
-  const auto d1     = refTraj.direction( zState );
-  const auto c1     = refTraj.curvature( zState );
-  const auto mat    = std::array{d0.mag2(), -d1.Dot( d0 ), d1.mag2() - dist.Dot( c1 )};
-  auto       mu     = std::array{-dist.Dot( d0 ), dist.Dot( d1 )};
-  const auto decomp = ROOT::Math::CholeskyDecomp<double, 2>{mat.data()};
+  // Determine the two points on the measurement and reference trajectory
+  // between which we minimize the distance between the trajectories
+  // ---> Yields Distance of Closest Approach (DOCA)
+  //
+  // This function is mainly used from the track fit.
+  // It provides extra output and handles potential local curvature of refTraj
+  // If you only need the doca or the location of the two points for two straight trajectories
+  // have a look what's defined in LHCb/Kernel/LHCbMath/GeomFun.h
+  //
+  // The curvature of the reference trajectory is evaluted at zState and treated as
+  // constant in the proximity of zState. This assumption only holds,
+  // if the final point on the reference trajectory is in fact close to zState.
+  // The final point is usually so close that the curvature term doesn't matter much
+  // which is why most (all?) configurations call this function with a refTraj that
+  // doesn't account for the local B-field vector
+  //
+  // The distance of closest approach (doca) vector is defined as (p_r - p_m)
+  // where p_r and p_m are points on the reference and measurement trajectory.
+  // The below answers: if I start at a point p on each trajectory,
+  // how far along their respective direction d do I have to go,
+  // to reach the points p_r or p_m respectively.
+  // The result to this called mu[0] and mu[1] below
+  // This then defines the two points as:
+  // p_r = p_r_zstate - mu[0] * d_r (+ potentially terms for curvature)
+  // p_m = p_m_0 - mu[1] * d_m
+
+  // this is the distance between our starting points p_r_zstate and p_m_0
+  auto const dist = refTraj.position( zState ) - m.trajectory().position( 0 );
+  auto const d_m  = m.trajectory().direction( 0 );
+  auto const d_r  = refTraj.direction( zState );
+  // extra term corresponding to curvature if the reference traj includes treatment of B-field
+  // normally we don't do that, thus c_r is just a vector of 0s
+  auto const c_r = refTraj.curvature( zState );
+  // now we solve the linear system for mu
+  auto const mat    = std::array{d_m.mag2(), -d_r.Dot( d_m ), d_r.mag2() - dist.Dot( c_r )};
+  auto       mu     = std::array{-dist.Dot( d_m ), dist.Dot( d_r )};
+  auto const decomp = ROOT::Math::CholeskyDecomp<double, 2>{mat.data()};
   if ( !decomp.Solve( mu ) ) throw std::runtime_error( "singular matrix" );
-  zState += mu[1];
-  dist += mu[0] * d0 - mu[1] * ( d1 + 0.5 * mu[1] * c1 );
-  // Set up the vector onto which we project everything. This should
-  // actually be parallel to dist.
-  auto unitPoca = d0.Cross( d1 + mu[1] * c1 ).Unit();
+
+  // we move to the correct point on the reference trajectory
+  auto const p_r = zState - mu[1];
+  // we started at 0 so m_r = 0 - mu[0]
+  auto const m_r = -mu[0];
+
+  // distance vector between the two points of closest approach
+  auto const doca_vec = refTraj.position( p_r ) - m.trajectory().position( m_r );
+
+  // vector onto which we project doca_vec to obtain a signed real value for the doca.
+  auto const unitPoca = d_m.Cross( refTraj.direction( p_r ) ).Unit();
+
   // compute the projection matrix from parameter space onto the (signed!) unit
-  Gaudi::TrackProjectionMatrix H = dual( unitPoca ) * refTraj.derivative( zState );
-  return {zState, mu[0], unitPoca.Dot( dist ), std::move( H ), std::move( unitPoca )};
+  Gaudi::TrackProjectionMatrix H = dual( unitPoca ) * refTraj.derivative( p_r );
+
+  return {p_r, m_r, unitPoca.Dot( doca_vec ), std::move( H ), std::move( unitPoca )};
 }

Tr/TrackTools/src/MeasurementProviderProjector.cpp

@@ -29,18 +29,6 @@ StatusCode MeasurementProviderProjector::initialize() {
   return extends::initialize().andThen( [&] { m_pIMF = svc<IMagneticFieldSvc>( "MagneticFieldSvc", true ); } );
 }
 
-// trivial helpers to make code clearer...
-namespace {
-  typedef Gaudi::Matrix1x3 DualVector;
-
-  DualVector dual( const Gaudi::XYZVector& v ) {
-    DualVector d;
-    v.GetCoordinates( d.Array() );
-    return d;
-  }
-
-} // namespace
-
 MeasurementProviderProjector::InternalProjectResult
 MeasurementProviderProjector::internal_project( const LHCb::StateVector& statevector,
                                                 const LHCb::Measurement& meas ) const {
@@ -102,33 +90,3 @@ StatusCode MeasurementProviderProjector::project( const LHCb::StateVector& state
   } catch ( StatusCode scr ) { sc = scr; }
   return sc;
 }
-
-//-----------------------------------------------------------------------------
-/// Derivatives wrt.the measurement's alignment...
-//-----------------------------------------------------------------------------
-MeasurementProviderProjector::Derivatives
-MeasurementProviderProjector::alignmentDerivatives( const LHCb::StateVector& statevector, const LHCb::Measurement& meas,
-                                                    const Gaudi::XYZPoint& pivot ) const {
-  auto       result = internal_project( statevector, meas );
-  DualVector unit   = dual( result.unitPocaVector );
-
-  // Calculate the derivative of the poca on measTraj to alignment parameters.
-  // Only non-zero elements:
-  Gaudi::XYZVector                  arm = meas.trajectory().position( result.sMeas ) - pivot;
-  ROOT::Math::SMatrix<double, 3, 6> dPosdAlpha;
-  // Derivative to translation
-  dPosdAlpha( 0, 0 ) = dPosdAlpha( 1, 1 ) = dPosdAlpha( 2, 2 ) = 1;
-  // Derivative to rotation around x-axis
-  dPosdAlpha( 1, 3 ) = -arm.z();
-  dPosdAlpha( 2, 3 ) = arm.y();
-  // Derivative to rotation around y-axis
-  dPosdAlpha( 0, 4 ) = arm.z();
-  dPosdAlpha( 2, 4 ) = -arm.x();
-  // Derivative to rotation around z-axis
-  dPosdAlpha( 0, 5 ) = -arm.y();
-  dPosdAlpha( 1, 5 ) = arm.x();
-
-  return unit * dPosdAlpha;
-  // compute the projection matrix from parameter space onto the (signed!) unit
-  // return unit*AlignTraj( meas.trajectory(), pivot ).derivative( m_sMeas );
-}

Tr/TrackTools/src/MeasurementProviderProjector.h

@@ -50,12 +50,6 @@ public:
   /// reset internal state of the provider, if any
   void reset() override{};
 
-  /// Retrieve the derivative of the residual wrt. the alignment parameters
-  /// of the measurement. The details of the alignment transformation are
-  /// defined in AlignTraj.
-  virtual Derivatives alignmentDerivatives( const LHCb::StateVector& state, const LHCb::Measurement& meas,
-                                            const Gaudi::XYZPoint& pivot ) const;
-
   StatusCode project( const LHCb::StateVector& statevector, LHCb::Measurement& meas, Gaudi::TrackProjectionMatrix& H,
                       double& residual, double& errMeasure, Gaudi::XYZVector& unitPocaVector, double& doca,
                       double& sMeas ) const;