new Thor functions for BandQ triggers

Phys/FunctorCore/include/Functors/Combination.h

@@ -51,6 +51,18 @@ namespace Functors::detail {
   private:
     std::optional<DistanceCalculator> m_dist_calc;
   };
+
+  //  template <typename DistanceCalculator, typename T, T N, T M>
+  template <typename T, T N, T M>
+  struct CosAngleBetweenDecayProducts : public Function {
+    static_assert( N >= 1 && M >= 1, "Indices start from 1 for LoKi compatibility." );
+    CosAngleBetweenDecayProducts( std::integral_constant<T, N>, std::integral_constant<T, M> ) {}
+
+    template <typename CombinationType>
+    auto operator()( CombinationType const& combination ) const {
+      return combination.template cos_angle_prod<N - 1, M - 1>();
+    }
+  };
 } // namespace Functors::detail
 
 /** @namespace Functors::Combination
@@ -70,6 +82,12 @@ namespace Functors::Combination {
     return detail::DistanceOfClosestApproachChi2<DistanceCalculator, T, N, M>( {}, {} );
   }
 
+  template <typename T, T N, T M>
+  auto CosAngleBetweenDecayProducts( std::integral_constant<T, N>, std::integral_constant<T, M> ) {
+    // This needs a wrapper function so that `T`, `N` and `M` can be deduced
+    return detail::CosAngleBetweenDecayProducts<T, N, M>( {}, {} );
+  }
+
   template <typename DistanceCalculator = detail::DefaultDistanceCalculator_t>
   struct MaxDistanceOfClosestApproach : public Function {
     void bind( TopLevelInfo& top_level ) { m_dist_calc.emplace( top_level.algorithm() ); }

Phys/FunctorCore/include/Functors/Composite.h

@@ -101,6 +101,42 @@ namespace Functors::detail {
     }
   };
 
+  /** BPVVDZ */
+  /** origin version in Loki defined in PHYS/PHYS_v25r1/Phys/LoKiPhys/src/Particles20.cpp,
+  VertexZDistanceWithTheBestPV function */
+  struct DeltaZToVertex : public Function {
+    /** This allows some error messages to be customised. It is not critical. */
+    static constexpr auto name() { return "DeltaZToVertex"; }
+
+    template <typename VContainer, typename Particle>
+    auto operator()( VContainer const& vertices, Particle const& composite ) const {
+      // Get the position of associated PV
+      auto const& bestPV_position = getBestPVPosition( composite, vertices );
+      // Get the position of end vertex
+      auto const& decay_vertex = composite.endVertex();
+
+      return ( decay_vertex->position() - bestPV_position ).z();
+      // return  decay_vertex->position().z() - bestPV_position.z();
+    }
+  };
+
+  /** BPVVDRHO */
+  /** origin version in Loki defined in PHYS/PHYS_v25r1/Phys/LoKiPhys/src/Particles20.cpp,
+  VertexRhoDistanceWithTheBestPV function */
+  struct DeltaRhoToVertex : public Function {
+    /** This allows some error messages to be customised. It is not critical. */
+    static constexpr auto name() { return "DeltaRhoToVertex"; }
+
+    template <typename VContainer, typename Particle>
+    auto operator()( VContainer const& vertices, Particle const& composite ) const {
+      // Get the position of associated PV
+      auto const& bestPV_position = getBestPVPosition( composite, vertices );
+      // Get the position of end vertex
+      auto const& decay_vertex = composite.endVertex();
+      return ( decay_vertex->position() - bestPV_position ).Rho();
+    }
+  };
+
   /** BPVDIRA */
   struct CosDirectionAngleToVertex : public Function {
     /** This allows some error messages to be customised. It is not critical. */
@@ -122,7 +158,9 @@ namespace Functors::detail {
       }();
 
       // The flight direction vector
-      auto const flight = Sel::Utils::endVertexPos( composite ) - getBestPVPosition( composite, vertices );
+      auto const  flight          = Sel::Utils::endVertexPos( composite ) - getBestPVPosition( composite, vertices );
+      auto const& bestPV_position = getBestPVPosition( composite, vertices );
+      auto const& decay_vertex_v2 = Sel::Utils::endVertexPos( composite );
       using std::sqrt;
       return Sel::vector::dot{}( mom, flight ) /
              sqrt( Sel::vector::magnitude2{}(flight)*Sel::vector::magnitude2{}( mom ) );
@@ -192,6 +230,20 @@ namespace Functors::detail {
     Functors::detail::DefaultLifetimeFitter_t m_lifetime_calc;
   };
 
+  struct ComputeDecayLengthSignificance : public Function {
+    /** This allows some error messages to be customised. It is not critical. */
+    static constexpr auto name() { return "ComputeDecayLengthSignificance"; }
+
+    template <typename VContainer, typename Composite>
+    auto operator()( VContainer const& vertices, Composite const& composite ) const {
+      auto const& bestPV = Sel::getBestPV( composite, vertices );
+      return m_lifetime_calc.DecayLengthSignificance( bestPV, composite );
+    }
+
+  private:
+    Functors::detail::DefaultLifetimeFitter_t m_lifetime_calc;
+  };
+
 } // namespace Functors::detail
 
 namespace Functors::Composite {
@@ -207,6 +259,17 @@ namespace Functors::Composite {
         std::move( vertex_location )};
   }
 
+  /** @brief Z component of flight distance of the given composite. */
+  template <typename VContainer = detail::DefaultPVContainer_t>
+  auto DeltaZToVertex( std::string vertex_location ) {
+    return detail::DataDepWrapper<Function, detail::DeltaZToVertex, VContainer>{std::move( vertex_location )};
+  }
+
+  /** @brief Rho component of flight distance of the given composite. */
+  template <typename VContainer = detail::DefaultPVContainer_t>
+  auto DeltaRhoToVertex( std::string vertex_location ) {
+    return detail::DataDepWrapper<Function, detail::DeltaRhoToVertex, VContainer>{std::move( vertex_location )};
+  }
   /** @brief Calculate the cosine of the angle between the momentum vector of
    *         the composite and the vector connecting its decay vertex to the
    *         "best" one of the given vertices.
@@ -319,4 +382,11 @@ namespace Functors::Composite {
   auto ComputeLifetime( std::string vertex_location ) {
     return detail::DataDepWrapper<Function, detail::ComputeLifetime, VContainer>{std::move( vertex_location )};
   }
+
+  /** @brief Calculate the lifetime of the particle with respect to the PV. */
+  template <typename VContainer = detail::DefaultPVContainer_t>
+  auto ComputeDecayLengthSignificance( std::string vertex_location ) {
+    return detail::DataDepWrapper<Function, detail::ComputeDecayLengthSignificance, VContainer>{
+        std::move( vertex_location )};
+  }
 } // namespace Functors::Composite

Phys/FunctorCore/python/Functors/__init__.py

@@ -149,6 +149,17 @@ DOCACHI2 = Functor(
     TemplateParams=[('DistanceCalculator',
                      'Distance calculator implementation to use.')],
     AllowMultiplePositionalArguments=True)
+ALV = Functor(
+    'ALV',
+    "Combination::CosAngleBetweenDecayProducts",
+    "Combination.h",
+    """Compute the cosine value of angle between two decay products.
+    """,
+    Params=[('Child1',
+             'Index [starting from 1] of the first child to consider.', int),
+            ('Child2',
+             'Index [starting from 1] of the second child to consider.', int)],
+    AllowMultiplePositionalArguments=True)
 MAXDOCA = Functor(
     'DOCA', "Combination::MaxDistanceOfClosestApproach", "Combination.h",
     "Compute the maximum pairwise distance of closest approach between members of a combination."
@@ -234,6 +245,26 @@ BPVFDCHI2 = Functor(
     Params=[('Vertices', 'TES location of input [primary] vertices.',
              DataHandle)],
     TemplateParams=[('VerticesType', 'Input vertex container type')])
+BPVVDZ = Functor(
+    'BPVVDZ',
+    'Composite::DeltaZToVertex',
+    'Composite.h',
+    '''Return the z component of flight distance w.r.t. the associated vertex, assuming
+    that it is from the container of vertices that is passed.
+    ''',
+    Params=[('Vertices', 'TES location of input [primary] vertices.',
+             DataHandle)],
+    TemplateParams=[('VerticesType', 'Input vertex container type')])
+BPVVDRHO = Functor(
+    'BPVVDRHO',
+    'Composite::DeltaRhoToVertex',
+    'Composite.h',
+    '''Return the rho component of flight distance w.r.t. the associated vertex, assuming
+    that it is from the container of vertices that is passed.
+    ''',
+    Params=[('Vertices', 'TES location of input [primary] vertices.',
+             DataHandle)],
+    TemplateParams=[('VerticesType', 'Input vertex container type')])
 BPVLTIME = Functor(
     'BPVLTIME',
     'Composite::ComputeLifetime',
@@ -246,6 +277,18 @@ BPVLTIME = Functor(
     TemplateParams=[
         ('VerticesType', 'Input vertex container type'),
     ])
+BPVDLS = Functor(
+    'BPVDLS',
+    'Composite::ComputeDecayLengthSignificance',
+    'Composite.h',
+    '''Return the decay length significance w.r.t. the associated vertex, assuming
+    that it is from the container of vertices that is passed. If no association
+    is available, compute one.''',
+    Params=[('Vertices', 'TES location of input [primary] vertices.',
+             DataHandle)],
+    TemplateParams=[
+        ('VerticesType', 'Input vertex container type'),
+    ])
 RUNNUMBER = Functor(
     'RUNNUMBER',
     'TES::RunNumber',

Phys/SelTools/include/SelTools/DistanceCalculator.h

@@ -222,6 +222,51 @@ namespace Sel {
     // LifetimeFitter( Gaudi::Algorithm* ) {}
   private:
   public:
+    /** @fn   DecayLengthSignificance
+     *  @brief Calculate the significance of a non-zero decay length.
+     *
+     */
+    template <typename Particle, typename VContainer>
+    auto DecayLengthSignificance( VContainer const& primary, Particle const& particle ) const {
+
+      // Calculate the distance between the particle and the vertex we hold.
+      Gaudi::XYZVector f = particle.endVertex()->position() - primary.position();
+      Gaudi::Vector3   flight( f.x(), f.y(), f.z() );
+
+      // Get the 3 momentum vectors for following calculations.
+      Gaudi::XYZVector p3    = particle.momentum().Vect();
+      double           p3mag = p3.R();
+      Gaudi::XYZVector d     = p3.Unit();
+      Gaudi::Vector3   dir( d.x(), d.y(), d.z() );
+
+      // Get the covariance matrix of the vertex.
+      Gaudi::SymMatrix3x3 W = primary.covMatrix();
+
+      // Get the particle covariance matrix.
+      const Gaudi::SymMatrix7x7 pCov = particle.covMatrix();
+
+      // Project the momentum of the particle onto its distance to the vertex
+      double a = ROOT::Math::Dot( dir, flight ) / p3mag;
+
+      // Update the covariance matrix
+      for ( size_t row = 0; row < 3; ++row ) {
+        for ( size_t col = 0; col <= row; ++col ) {
+          W( row, col ) +=
+              pCov( row, col ) + a * a * pCov( row + 3, col + 3 ) - a * ( pCov( row + 3, col ) + pCov( col + 3, row ) );
+        }
+      }
+      // int OK = W.Invert();
+      //      if( !OK ) {
+      //          Warning( "Error inverting the updated covariance matrix." ).ignore();
+      //      }
+      // Calculate the required numbers.
+      double halfdChi2dLam2 = ROOT::Math::Similarity( W, dir );
+      double decayLength    = ROOT::Math::Dot( dir, W * flight ) / halfdChi2dLam2;
+      double decayLengthErr = std::sqrt( 1 / halfdChi2dLam2 );
+
+      return decayLength / decayLengthErr;
+    }
+
     /** @fn    Lifetime
      *  @brief Calculate the lifetime between the particle and the best PV.
      *