Complete RecoDecay. Update related tables and tasks. (#4256)

Implement calculation of decay length uncertainty.
Replace DCA with Chi2.
Rename tasks and move to PWGHF.
Remove obsolete files.
Minor fixes.

Author: vkucera

Analysis/DataModel/CMakeLists.txt

@@ -9,8 +9,7 @@
 # submit itself to any jurisdiction.
 
 o2_add_library(AnalysisDataModel
-               SOURCES src/dummy.cxx
-                       src/ParamBase.cxx
+               SOURCES src/ParamBase.cxx
                        src/DetectorResponse.cxx
                        src/PIDTOF.cxx
                        src/PIDTPC.cxx

Analysis/DataModel/include/Analysis/RecoDecay.h

@@ -56,7 +56,7 @@ class RecoDecay
   /// \param num  a number of arbitrary type
   /// \return number squared
   template <typename T>
-  auto static sq(T num)
+  static auto sq(T num)
   {
     return (double)num * (double)num;
   }
@@ -128,16 +128,16 @@ class RecoDecay
   /// \param point1,point2  {x, y, z} coordinates of points
   /// \return 3D distance between two points
   template <typename T, typename U>
-  static auto distance(const array<T, 3>& point1, const array<U, 3>& point2)
+  static auto distance(const T& point1, const U& point2)
   {
     return sqrtSumOfSquares(point1[0] - point2[0], point1[1] - point2[1], point1[2] - point2[2]);
   }
 
   /// Calculates 2D {x, y} distance between two points.
   /// \param point1,point2  {x, y, z} or {x, y} coordinates of points
   /// \return 2D {x, y} distance between two points
-  template <std::size_t N, std::size_t O, typename T, typename U>
-  static auto distanceXY(const array<T, N>& point1, const array<U, O>& point2)
+  template <typename T, typename U>
+  static auto distanceXY(const T& point1, const U& point2)
   {
     return sqrtSumOfSquares(point1[0] - point2[0], point1[1] - point2[1]);
   }
@@ -279,8 +279,8 @@ class RecoDecay
   /// \note AliRoot version
   /// \param arrMom  array of two 3-momentum arrays
   /// \param arrMass  array of two masses (in the same order as arrMom)
-  /// \param iProng  index of the prong
   /// \param mTot  assumed mass of mother particle
+  /// \param iProng  index of the prong
   /// \return cosine of θ* of the i-th prong under the assumption of the invariant mass
   template <typename T, typename U, typename V>
   static auto CosThetaStarA(const array<array<T, 3>, 2>& arrMom, const array<U, 2>& arrMass, V mTot, int iProng)
@@ -294,7 +294,7 @@ class RecoDecay
     // p* = √[(M^2 - m1^2 - m2^2)^2 - 4 m1^2 m2^2]/2M
     // Lorentz transformation of the longitudinal momentum of the prong into the detector frame:
     // p_L,i = γ (p*_L,i + β E*_i)
-    // p*_L,i = p_L,i/γ - β E*_i)
+    // p*_L,i = p_L,i/γ - β E*_i
     // cos(θ*_i) = (p_L,i/γ - β E*_i)/p*
     return (dotProd(arrMom[iProng], pVecTot) / (pTot * gamma) - beta * E(pStar, arrMass[iProng])) / pStar;
   }

Analysis/DataModel/include/Analysis/SecondaryVertex.h

@@ -56,6 +56,7 @@ DECLARE_SOA_COLUMN(XSecondaryVertex, xSecondaryVertex, float);
 DECLARE_SOA_COLUMN(YSecondaryVertex, ySecondaryVertex, float);
 DECLARE_SOA_COLUMN(ZSecondaryVertex, zSecondaryVertex, float);
 DECLARE_SOA_DYNAMIC_COLUMN(RSecondaryVertex, rSecondaryVertex, [](float xVtxS, float yVtxS) { return RecoDecay::sqrtSumOfSquares(xVtxS, yVtxS); });
+DECLARE_SOA_COLUMN(Chi2PCA, chi2PCA, float); // sum of (non-weighted) distances of the secondary vertex to its prongs
 // prong properties
 DECLARE_SOA_COLUMN(PxProng0, pxProng0, float);
 DECLARE_SOA_COLUMN(PyProng0, pyProng0, float);
@@ -108,7 +109,6 @@ DECLARE_SOA_DYNAMIC_COLUMN(Ct, ct, [](float xVtxP, float yVtxP, float zVtxP, flo
 // specific 2-prong decay properties
 namespace hf_cand_prong2
 {
-DECLARE_SOA_COLUMN(DCA, dca, float);
 DECLARE_SOA_EXPRESSION_COLUMN(Px, px, float, 1.f * aod::hf_cand::pxProng0 + 1.f * aod::hf_cand::pxProng1);
 DECLARE_SOA_EXPRESSION_COLUMN(Py, py, float, 1.f * aod::hf_cand::pyProng0 + 1.f * aod::hf_cand::pyProng1);
 DECLARE_SOA_EXPRESSION_COLUMN(Pz, pz, float, 1.f * aod::hf_cand::pzProng0 + 1.f * aod::hf_cand::pzProng1);
@@ -122,43 +122,43 @@ DECLARE_SOA_DYNAMIC_COLUMN(CosThetaStar, cosThetaStar, [](float px0, float py0,
 // D0(bar) → π± K∓
 
 template <typename T>
-auto CtD0(T candidate)
+auto CtD0(const T& candidate)
 {
   return candidate.ct(RecoDecay::getMassPDG(421));
 }
 
 template <typename T>
-auto YD0(T candidate)
+auto YD0(const T& candidate)
 {
   return candidate.y(RecoDecay::getMassPDG(421));
 }
 
 template <typename T>
-auto ED0(T candidate)
+auto ED0(const T& candidate)
 {
   return candidate.e(RecoDecay::getMassPDG(421));
 }
 
 template <typename T>
-auto InvMassD0(T candidate)
+auto InvMassD0(const T& candidate)
 {
   return candidate.m(array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kKPlus)});
 }
 
 template <typename T>
-auto InvMassD0bar(T candidate)
+auto InvMassD0bar(const T& candidate)
 {
   return candidate.m(array{RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kPiPlus)});
 }
 
 template <typename T>
-auto CosThetaStarD0(T candidate)
+auto CosThetaStarD0(const T& candidate)
 {
   return candidate.cosThetaStar(array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kKPlus)}, RecoDecay::getMassPDG(421), 1);
 }
 
 template <typename T>
-auto CosThetaStarD0bar(T candidate)
+auto CosThetaStarD0bar(const T& candidate)
 {
   return candidate.cosThetaStar(array{RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kPiPlus)}, RecoDecay::getMassPDG(421), 0);
 }
@@ -168,6 +168,7 @@ DECLARE_SOA_TABLE(HfCandBase, "AOD", "HFCANDBASE",
                   collision::PosX, collision::PosY, collision::PosZ,
                   hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex,
                   hf_cand::ErrorDecayLength, hf_cand::ErrorDecayLengthXY,
+                  hf_cand::Chi2PCA,
                   /* dynamic columns */
                   hf_cand::RSecondaryVertex<hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex>,
                   hf_cand::ImpactParameterNormalised0<hf_cand::ImpactParameter0, hf_cand::ErrorImpactParameter0>,
@@ -188,7 +189,6 @@ DECLARE_SOA_TABLE(HfCandBase, "AOD", "HFCANDBASE",
                   hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1,
                   hf_cand::ImpactParameter0, hf_cand::ImpactParameter1,
                   hf_cand::ErrorImpactParameter0, hf_cand::ErrorImpactParameter1,
-                  hf_cand_prong2::DCA,
                   /* dynamic columns */
                   hf_cand_prong2::ImpactParameterProduct<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1>,
                   hf_cand_prong2::M<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1>,

Analysis/DataModel/include/Analysis/SecondaryVertexHF.h

@@ -7,21 +7,25 @@
 // In applying this license CERN does not waive the privileges and immunities
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
-#include "ReconstructionDataFormats/Track.h"
-#include "ReconstructionDataFormats/Vertex.h"
 
 /// \file trackUtilities.h
-/// \brief Utilities for manipulating track parameters
+/// \brief Utilities for manipulating parameters of tracks and vertices
 ///
 /// \author Vít Kučera <vit.kucera@cern.ch>, CERN
 
 #ifndef O2_ANALYSIS_TRACKUTILITIES_H_
 #define O2_ANALYSIS_TRACKUTILITIES_H_
 
+#include <cmath>
+
+#include "ReconstructionDataFormats/Track.h"
+#include "ReconstructionDataFormats/Vertex.h"
+#include "Analysis/RecoDecay.h"
+
 /// Extracts track parameters and covariance matrix from a track.
 /// \return o2::track::TrackParCov
 template <typename T>
-auto getTrackParCov(T track)
+auto getTrackParCov(const T& track)
 {
   std::array<float, 5> arraypar = {track.y(), track.z(), track.snp(),
                                    track.tgl(), track.signed1Pt()};
@@ -37,11 +41,47 @@ auto getTrackParCov(T track)
 /// Extracts primary vertex position and covariance matrix from a collision.
 /// \return o2::dataformats::VertexBase
 template <typename T>
-auto getPrimaryVertex(T collision)
+auto getPrimaryVertex(const T& collision)
 {
   Point3D<float> vtxXYZ(collision.posX(), collision.posY(), collision.posZ());
   std::array<float, 6> vtxCov{collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()};
   return o2::dataformats::VertexBase{std::move(vtxXYZ), std::move(vtxCov)};
 }
 
+/// Calculates direction of one point w.r.t another point.
+/// \param point1,point2  points with {x, y, z} coordinates accessible by index
+/// \param phi  azimuth angle in the {x, y} plane (taken w.r.t. the x-axis towards the y-axis)
+/// \param theta  angle w.r.t the {x, y} plane towards the z-axis
+/// \return phi,theta of point2 w.r.t. point1
+template <typename T, typename U, typename V, typename W>
+void getPointDirection(const T& point1, const U& point2, V& phi, W& theta)
+{
+  phi = std::atan2(point2[1] - point1[1], point2[0] - point1[0]);
+  //auto x1 = point1[0] * std::cos(phi) + point1[1] * std::sin(phi);
+  //auto x2 = point2[0] * std::cos(phi) + point2[1] * std::sin(phi);
+  //theta = std::atan2(point2[2] - point1[2], x2 - x1);
+  theta = std::atan2(point2[2] - point1[2], RecoDecay::distanceXY(point1, point2));
+}
+
+/// Calculates the XX element of a XYZ covariance matrix after rotation of the coordinate system
+/// by phi around the z-axis and by minus theta around the new y-axis.
+/// \param matrix  matrix
+/// \param phi  azimuth angle in the {x, y} plane (taken w.r.t. the x-axis towards the y-axis)
+/// \param theta  angle w.r.t the {x, y} plane towards the z-axis
+/// \return XX element of the rotated covariance matrix
+template <typename T, typename U, typename V>
+auto getRotatedCovMatrixXX(const T& matrix, U phi, V theta)
+{
+  auto cp = std::cos(phi);
+  auto sp = std::sin(phi);
+  auto ct = std::cos(theta);
+  auto st = std::sin(theta);
+  return matrix[0] * cp * cp * ct * ct        // covXX
+         + matrix[1] * 2. * cp * sp * ct * ct // covXY
+         + matrix[2] * sp * sp * ct * ct      // covYY
+         + matrix[3] * 2. * cp * ct * st      // covXZ
+         + matrix[4] * 2. * sp * ct * st      // covYZ
+         + matrix[5] * st * st;               // covZZ
+}
+
 #endif // O2_ANALYSIS_TRACKUTILITIES_H_