diff --git a/Analysis/Core/include/Analysis/RecoDecay.h b/Analysis/Core/include/Analysis/RecoDecay.h
index 530b920724d08..87ef939b8a25c 100644
--- a/Analysis/Core/include/Analysis/RecoDecay.h
+++ b/Analysis/Core/include/Analysis/RecoDecay.h
@@ -101,7 +101,7 @@ class RecoDecay
   }
 
   /// Calculates scalar product of vectors.
-  /// \note Promotes numbers to double before squaring to avoid precision loss in float multiplication.
+  /// \note Promotes numbers to double to avoid precision loss in float multiplication.
   /// \param N  dimension
   /// \param vec1,vec2  vectors
   /// \return scalar product
@@ -115,6 +115,19 @@ class RecoDecay
     return res;
   }
 
+  /// FIXME: probably cross and dot products should be in some utility class
+  /// Calculates cross product of vectors in three dimensions.
+  /// \note Promotes numbers to double to avoid precision loss in float multiplication.
+  /// \param vec1,vec2  vectors
+  /// \return cross-product vector
+  template <typename T, typename U>
+  static array<double, 3> crossProd(const array<T, 3>& vec1, const array<U, 3>& vec2)
+  {
+    return array<double, 3>{((double)vec1[1] * (double)vec2[2]) - ((double)vec1[2] * (double)vec2[1]),
+                            ((double)vec1[2] * (double)vec2[0]) - ((double)vec1[0] * (double)vec2[2]),
+                            ((double)vec1[0] * (double)vec2[1]) - ((double)vec1[1] * (double)vec2[0])};
+  }
+
   /// Calculates magnitude squared of a vector.
   /// \param N  dimension
   /// \param vec  vector
@@ -388,6 +401,76 @@ class RecoDecay
     return std::sqrt(M2(args...));
   }
 
+  // Calculation of topological quantities
+
+  /// Calculates impact parameter in the bending plane of the particle w.r.t. a point
+  /// \param point  {x, y, z} position of the point
+  /// \param posSV  {x, y, z} position of the secondary vertex
+  /// \param mom  {x, y, z} particle momentum array
+  /// \return impact parameter in {x, y}
+  template <typename T, typename U, typename V>
+  static double ImpParXY(const T& point, const U& posSV, const array<V, 3>& mom)
+  {
+    // Ported from AliAODRecoDecay::ImpParXY
+    auto flightLineXY = array{posSV[0] - point[0], posSV[1] - point[1]};
+    auto k = dotProd(flightLineXY, array{mom[0], mom[1]}) / Pt2(mom);
+    auto dx = flightLineXY[0] - k * (double)mom[0];
+    auto dy = flightLineXY[1] - k * (double)mom[1];
+    auto absImpPar = sqrtSumOfSquares(dx, dy);
+    auto flightLine = array{posSV[0] - point[0], posSV[1] - point[1], posSV[2] - point[2]};
+    auto cross = crossProd(mom, flightLine);
+    return (cross[2] > 0. ? absImpPar : -1. * absImpPar);
+  }
+
+  /// Calculates the difference between measured and expected track impact parameter
+  /// normalized to its uncertainty
+  /// \param decLenXY decay lenght in {x, y} plane
+  /// \param errDecLenXY error on decay lenght in {x, y} plane
+  /// \param momMother {x, y, z} or {x, y} candidate momentum array
+  /// \param impParProng prong impact parameter
+  /// \param errImpParProng error on prong impact parameter
+  /// \param momProng {x, y, z} or {x, y} prong momentum array
+  /// \return normalized difference between expected and observed impact parameter
+  template <std::size_t N, std::size_t M, typename T, typename U, typename V, typename W, typename X, typename Y>
+  static double normImpParMeasMinusExpProng(T decLenXY, U errDecLenXY, const array<V, N>& momMother, W impParProng,
+                                            X errImpParProng, const array<Y, M>& momProng)
+  {
+    // Ported from AliAODRecoDecayHF::Getd0MeasMinusExpProng adding normalization directly in the function
+    auto sinThetaP = ((double)momProng[0] * (double)momMother[1] - (double)momProng[1] * (double)momMother[0]) /
+                     (Pt(momProng) * Pt(momMother));
+    auto diff = impParProng - (double)decLenXY * sinThetaP;
+    auto errImpParExpProng = (double)errDecLenXY * sinThetaP;
+    auto errDiff = sqrtSumOfSquares(errImpParProng, errImpParExpProng);
+    return (errDiff > 0. ? diff / errDiff : 0.);
+  }
+
+  /// Calculates maximum normalized difference between measured and expected impact parameter of candidate prongs
+  /// \param posPV {x, y, z} or {x, y} position of primary vertex
+  /// \param posSV {x, y, z} or {x, y} position of secondary vertex
+  /// \param errDecLenXY error on decay lenght in {x, y} plane
+  /// \param momMother {x, y, z} or {x, y} candidate momentum array
+  /// \param arrImpPar array of prong impact parameters
+  /// \param arrErrImpPar array of errors on prong impact parameter (same order as arrImpPar)
+  /// \param momMom array of {x, y, z} or {x, y} prong momenta (same order as arrImpPar)
+  /// \return maximum normalized difference between expected and observed impact parameters
+  template <std::size_t N, std::size_t M, std::size_t K, typename T, typename U, typename V, typename W, typename X,
+            typename Y, typename Z>
+  static double maxNormalisedDeltaIP(const T& posPV, const U& posSV, V errDecLenXY, const array<W, M>& momMother,
+                                     const array<X, N>& arrImpPar, const array<Y, N>& arrErrImpPar,
+                                     const array<array<Z, K>, N>& arrMom)
+  {
+    auto decLenXY = distanceXY(posPV, posSV);
+    double maxNormDeltaIP{0.};
+    for (auto iProng = 0; iProng < N; ++iProng) {
+      auto prongNormDeltaIP = normImpParMeasMinusExpProng(decLenXY, errDecLenXY, momMother, arrImpPar[iProng],
+                                                          arrErrImpPar[iProng], arrMom[iProng]);
+      if (std::abs(prongNormDeltaIP) > std::abs(maxNormDeltaIP)) {
+        maxNormDeltaIP = prongNormDeltaIP;
+      }
+    }
+    return maxNormDeltaIP;
+  }
+
   /// Returns particle mass based on PDG code.
   /// \param pdg  PDG code
   /// \return particle mass
diff --git a/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h b/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h
index 8ee186e65b55d..88049ea10ea9f 100644
--- a/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h
+++ b/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h
@@ -119,6 +119,7 @@ DECLARE_SOA_COLUMN(ErrorDecayLengthXY, errorDecayLengthXY, float);
 DECLARE_SOA_DYNAMIC_COLUMN(CPA, cpa, [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz) { return RecoDecay::CPA(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}, array{px, py, pz}); });
 DECLARE_SOA_DYNAMIC_COLUMN(CPAXY, cpaXY, [](float xVtxP, float yVtxP, float xVtxS, float yVtxS, float px, float py) { return RecoDecay::CPAXY(array{xVtxP, yVtxP}, array{xVtxS, yVtxS}, array{px, py}); });
 DECLARE_SOA_DYNAMIC_COLUMN(Ct, ct, [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz, double m) { return RecoDecay::Ct(array{px, py, pz}, RecoDecay::distance(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}), m); });
+DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterXY, impactParameterXY, [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz) { return RecoDecay::ImpParXY(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}, array{px, py, pz}); });
 } // namespace hf_cand
 
 // specific 2-prong decay properties
@@ -131,6 +132,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProduct, impactParameterProduct, [](fl
 DECLARE_SOA_DYNAMIC_COLUMN(M, m, [](float px0, float py0, float pz0, float px1, float py1, float pz1, const array<double, 2>& m) { return RecoDecay::M(array{array{px0, py0, pz0}, array{px1, py1, pz1}}, m); });
 DECLARE_SOA_DYNAMIC_COLUMN(M2, m2, [](float px0, float py0, float pz0, float px1, float py1, float pz1, const array<double, 2>& m) { return RecoDecay::M2(array{array{px0, py0, pz0}, array{px1, py1, pz1}}, m); });
 DECLARE_SOA_DYNAMIC_COLUMN(CosThetaStar, cosThetaStar, [](float px0, float py0, float pz0, float px1, float py1, float pz1, const array<double, 2>& m, double mTot, int iProng) { return RecoDecay::CosThetaStar(array{array{px0, py0, pz0}, array{px1, py1, pz1}}, m, mTot, iProng); });
+DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProngSqSum, impactParameterProngSqSum, [](float impParProng0, float impParProng1) { return RecoDecay::sumOfSquares(impParProng0, impParProng1); });
+DECLARE_SOA_DYNAMIC_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, [](float xVtxP, float yVtxP, float xVtxS, float yVtxS, float errDlxy, float pxM, float pyM, float ip0, float errIp0, float ip1, float errIp1, float px0, float py0, float px1, float py1) { return RecoDecay::maxNormalisedDeltaIP(array{xVtxP, yVtxP}, array{xVtxS, yVtxS}, errDlxy, array{pxM, pyM}, array{ip0, ip1}, array{errIp0, errIp1}, array{array{px0, py0}, array{px1, py1}}); });
 // MC matching result:
 // - bit 0: D0(bar) → π± K∓
 DECLARE_SOA_COLUMN(FlagMCMatchRec, flagMCMatchRec, uint8_t); // reconstruction level
@@ -218,6 +221,7 @@ DECLARE_SOA_TABLE(HfCandProng2Base, "AOD", "HFCANDP2BASE",
                   hf_cand_prong2::M2<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1>,
                   hf_cand_prong2::ImpactParameterProduct<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1>,
                   hf_cand_prong2::CosThetaStar<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1>,
+                  hf_cand_prong2::ImpactParameterProngSqSum<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1>,
                   /* dynamic columns that use candidate momentum components */
                   hf_cand::Pt<hf_cand_prong2::Px, hf_cand_prong2::Py>,
                   hf_cand::Pt2<hf_cand_prong2::Px, hf_cand_prong2::Py>,
@@ -227,6 +231,8 @@ DECLARE_SOA_TABLE(HfCandProng2Base, "AOD", "HFCANDP2BASE",
                   hf_cand::CPA<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
                   hf_cand::CPAXY<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand_prong2::Px, hf_cand_prong2::Py>,
                   hf_cand::Ct<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
+                  hf_cand::ImpactParameterXY<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
+                  hf_cand_prong2::MaxNormalisedDeltaIP<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ErrorDecayLengthXY, hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand::ImpactParameter0, hf_cand::ErrorImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ErrorImpactParameter1, hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PxProng1, hf_cand::PyProng1>,
                   hf_cand::Eta<hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
                   hf_cand::Y<hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
                   hf_cand::E<hf_cand_prong2::Px, hf_cand_prong2::Py, hf_cand_prong2::Pz>,
@@ -254,6 +260,8 @@ DECLARE_SOA_EXPRESSION_COLUMN(Py, py, float, 1.f * aod::hf_cand::pyProng0 + 1.f
 DECLARE_SOA_EXPRESSION_COLUMN(Pz, pz, float, 1.f * aod::hf_cand::pzProng0 + 1.f * aod::hf_cand::pzProng1 + 1.f * aod::hf_cand::pzProng2);
 DECLARE_SOA_DYNAMIC_COLUMN(M, m, [](float px0, float py0, float pz0, float px1, float py1, float pz1, float px2, float py2, float pz2, const array<double, 3>& m) { return RecoDecay::M(array{array{px0, py0, pz0}, array{px1, py1, pz1}, array{px2, py2, pz2}}, m); });
 DECLARE_SOA_DYNAMIC_COLUMN(M2, m2, [](float px0, float py0, float pz0, float px1, float py1, float pz1, float px2, float py2, float pz2, const array<double, 3>& m) { return RecoDecay::M2(array{array{px0, py0, pz0}, array{px1, py1, pz1}, array{px2, py2, pz2}}, m); });
+DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterProngSqSum, impactParameterProngSqSum, [](float impParProng0, float impParProng1, float impParProng2) { return RecoDecay::sumOfSquares(impParProng0, impParProng1, impParProng2); });
+DECLARE_SOA_DYNAMIC_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, [](float xVtxP, float yVtxP, float xVtxS, float yVtxS, float errDlxy, float pxM, float pyM, float ip0, float errIp0, float ip1, float errIp1, float ip2, float errIp2, float px0, float py0, float px1, float py1, float px2, float py2) { return RecoDecay::maxNormalisedDeltaIP(array{xVtxP, yVtxP}, array{xVtxS, yVtxS}, errDlxy, array{pxM, pyM}, array{ip0, ip1, ip2}, array{errIp0, errIp1, errIp2}, array{array{px0, py0}, array{px1, py1}, array{px2, py2}}); });
 // MC matching result:
 // - bit 0: D± → π± K∓ π±
 // - bit 1: Lc± → p± K∓ π±
@@ -329,6 +337,7 @@ DECLARE_SOA_TABLE(HfCandProng3Base, "AOD", "HFCANDP3BASE",
                   /* dynamic columns */
                   hf_cand_prong3::M<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, hf_cand::PxProng2, hf_cand::PyProng2, hf_cand::PzProng2>,
                   hf_cand_prong3::M2<hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, hf_cand::PxProng2, hf_cand::PyProng2, hf_cand::PzProng2>,
+                  hf_cand_prong3::ImpactParameterProngSqSum<hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ImpactParameter2>,
                   /* prong 2 */
                   hf_cand::ImpactParameterNormalised2<hf_cand::ImpactParameter2, hf_cand::ErrorImpactParameter2>,
                   hf_cand::PtProng2<hf_cand::PxProng2, hf_cand::PyProng2>,
@@ -343,6 +352,8 @@ DECLARE_SOA_TABLE(HfCandProng3Base, "AOD", "HFCANDP3BASE",
                   hf_cand::CPA<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
                   hf_cand::CPAXY<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand_prong3::Px, hf_cand_prong3::Py>,
                   hf_cand::Ct<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
+                  hf_cand::ImpactParameterXY<collision::PosX, collision::PosY, collision::PosZ, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ZSecondaryVertex, hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
+                  hf_cand_prong3::MaxNormalisedDeltaIP<collision::PosX, collision::PosY, hf_cand::XSecondaryVertex, hf_cand::YSecondaryVertex, hf_cand::ErrorDecayLengthXY, hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand::ImpactParameter0, hf_cand::ErrorImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ErrorImpactParameter1, hf_cand::ImpactParameter2, hf_cand::ErrorImpactParameter2, hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PxProng2, hf_cand::PyProng2>,
                   hf_cand::Eta<hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
                   hf_cand::Y<hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
                   hf_cand::E<hf_cand_prong3::Px, hf_cand_prong3::Py, hf_cand_prong3::Pz>,
diff --git a/Analysis/Tasks/PWGHF/taskDPlus.cxx b/Analysis/Tasks/PWGHF/taskDPlus.cxx
index e422476a556f0..cf2af649c6c49 100644
--- a/Analysis/Tasks/PWGHF/taskDPlus.cxx
+++ b/Analysis/Tasks/PWGHF/taskDPlus.cxx
@@ -16,6 +16,7 @@
 
 #include "Framework/runDataProcessing.h"
 #include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
 #include "Analysis/HFSecondaryVertex.h"
 #include "Analysis/HFCandidateSelectionTables.h"
 
@@ -26,23 +27,32 @@ using namespace o2::aod::hf_cand_prong3;
 
 /// D± analysis task
 struct TaskDPlus {
-  OutputObj<TH1F> hmass{TH1F("hmass", "3-prong candidates;inv. mass (#pi K #pi) (GeV/#it{c}^{2});entries", 500, 1.6, 2.1)};
-  OutputObj<TH1F> hptcand{TH1F("hptcand", "3-prong candidates;candidate #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)};
-  OutputObj<TH1F> hptprong0{TH1F("hptprong0", "3-prong candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)};
-  OutputObj<TH1F> hptprong1{TH1F("hptprong1", "3-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)};
-  OutputObj<TH1F> hptprong2{TH1F("hptprong2", "3-prong candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)};
-  OutputObj<TH1F> hdeclength{TH1F("declength", "3-prong candidates;decay length (cm);entries", 200, 0., 2.)};
-  OutputObj<TH1F> hdeclengthxy{TH1F("declengthxy", "3-prong candidates;decay length xy (cm);entries", 200, 0., 2.)};
-  OutputObj<TH1F> hd0Prong0{TH1F("hd0Prong0", "3-prong candidates;prong 0 DCAxy to prim. vertex (cm);entries", 100, -1., 1.)};
-  OutputObj<TH1F> hd0Prong1{TH1F("hd0Prong1", "3-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", 100, -1., 1.)};
-  OutputObj<TH1F> hd0Prong2{TH1F("hd0Prong2", "3-prong candidates;prong 2 DCAxy to prim. vertex (cm);entries", 100, -1., 1.)};
-  OutputObj<TH1F> hCt{TH1F("hCt", "3-prong candidates;proper lifetime (D^{#pm}) * #it{c} (cm);entries", 120, -20., 100.)};
-  OutputObj<TH1F> hCPA{TH1F("hCPA", "3-prong candidates;cosine of pointing angle;entries", 110, -1.1, 1.1)};
-  OutputObj<TH1F> hEta{TH1F("hEta", "3-prong candidates;candidate #it{#eta};entries", 100, -2., 2.)};
-  //OutputObj<TH1F> hselectionstatus{TH1F("selectionstatus", "3-prong candidates;selection status;entries", 5, -0.5, 4.5)};
-  OutputObj<TH1F> hImpParErr{TH1F("hImpParErr", "3-prong candidates;impact parameter error (cm);entries", 100, -1., 1.)};
-  OutputObj<TH1F> hDecLenErr{TH1F("hDecLenErr", "3-prong candidates;decay length error (cm);entries", 100, 0., 1.)};
-  OutputObj<TH1F> hDecLenXYErr{TH1F("hDecLenXYErr", "3-prong candidates;decay length xy error (cm);entries", 100, 0., 1.)};
+  HistogramRegistry registry{
+    "registry",
+    {
+      {"hMass", "3-prong candidates;inv. mass (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{350, 1.7, 2.05}}}},
+      {"hPt", "3-prong candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}},
+      {"hEta", "3-prong candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}},
+      {"hCt", "3-prong candidates;proper lifetime (D^{#pm}) * #it{c} (cm);entries", {HistType::kTH1F, {{120, -20., 100.}}}},
+      {"hDecayLength", "3-prong candidates;decay length (cm);entries", {HistType::kTH1F, {{200, 0., 2.}}}},
+      {"hDecayLengthXY", "3-prong candidates;decay length xy (cm);entries", {HistType::kTH1F, {{200, 0., 2.}}}},
+      {"hNormalisedDecayLengthXY", "3-prong candidates;norm. decay length xy;entries", {HistType::kTH1F, {{80, 0., 80.}}}},
+      {"hCPA", "3-prong candidates;cos. pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}},
+      {"hCPAxy", "3-prong candidates;cos. pointing angle xy;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}},
+      {"hImpactParameterXY", "3-prong candidates;impact parameter xy (cm);entries", {HistType::kTH1F, {{200, -1., 1.}}}},
+      {"hMaxNormalisedDeltaIP", "3-prong candidates;norm. IP;entries", {HistType::kTH1F, {{200, -20., 20.}}}},
+      {"hImpactParameterProngSqSum", "3-prong candidates;squared sum of prong imp. par. (cm^{2});entries", {HistType::kTH1F, {{100, 0., 1.}}}},
+      {"hDecayLengthError", "3-prong candidates;decay length error (cm);entries", {HistType::kTH1F, {{100, 0., 1.}}}},
+      {"hDecayLengthXYError", "3-prong candidates;decay length xy error (cm);entries", {HistType::kTH1F, {{100, 0., 1.}}}},
+      {"hImpactParameterError", "3-prong candidates;impact parameter error (cm);entries", {HistType::kTH1F, {{100, 0., 1.}}}},
+      {"hPtProng0", "3-prong candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}},
+      {"hPtProng1", "3-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}},
+      {"hPtProng2", "3-prong candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}},
+      {"hd0Prong0", "3-prong candidates;prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{100, -1., 1.}}}},
+      {"hd0Prong1", "3-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{100, -1., 1.}}}},
+      {"hd0Prong2", "3-prong candidates;prong 2 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{100, -1., 1.}}}}
+      //{"hSelectionStatus", "3-prong candidates;selection status;entries",  {HistType::kTH1F, {{5, -0.5, 4.5}}}}
+    }};
 
   Configurable<int> d_selectionFlagDPlus{"d_selectionFlagDPlus", 1, "Selection Flag for DPlus"};
 
@@ -52,25 +62,30 @@ struct TaskDPlus {
   void process(aod::HfCandProng3 const& candidates)
   {
     for (auto& candidate : candidates) {
-      hmass->Fill(InvMassDPlus(candidate));
-      hptcand->Fill(candidate.pt());
-      hptprong0->Fill(candidate.ptProng0());
-      hptprong1->Fill(candidate.ptProng1());
-      hptprong2->Fill(candidate.ptProng2());
-      hdeclength->Fill(candidate.decayLength());
-      hdeclengthxy->Fill(candidate.decayLengthXY());
-      hd0Prong0->Fill(candidate.impactParameter0());
-      hd0Prong1->Fill(candidate.impactParameter1());
-      hd0Prong2->Fill(candidate.impactParameter2());
-      hCt->Fill(CtDPlus(candidate));
-      hCPA->Fill(candidate.cpa());
-      hEta->Fill(candidate.eta());
-      //hselectionstatus->Fill(candidate.isSelDPlus());
-      hImpParErr->Fill(candidate.errorImpactParameter0());
-      hImpParErr->Fill(candidate.errorImpactParameter1());
-      hImpParErr->Fill(candidate.errorImpactParameter2());
-      hDecLenErr->Fill(candidate.errorDecayLength());
-      hDecLenXYErr->Fill(candidate.errorDecayLengthXY());
+      registry.get<TH1>("hMass")->Fill(InvMassDPlus(candidate));
+      registry.get<TH1>("hPt")->Fill(candidate.pt());
+      registry.get<TH1>("hEta")->Fill(candidate.eta());
+      registry.get<TH1>("hCt")->Fill(CtDPlus(candidate));
+      registry.get<TH1>("hDecayLength")->Fill(candidate.decayLength());
+      registry.get<TH1>("hDecayLengthXY")->Fill(candidate.decayLengthXY());
+      registry.get<TH1>("hNormalisedDecayLengthXY")->Fill(candidate.decayLengthXYNormalised());
+      registry.get<TH1>("hCPA")->Fill(candidate.cpa());
+      registry.get<TH1>("hCPAxy")->Fill(candidate.cpaXY());
+      registry.get<TH1>("hImpactParameterXY")->Fill(candidate.impactParameterXY());
+      registry.get<TH1>("hMaxNormalisedDeltaIP")->Fill(candidate.maxNormalisedDeltaIP());
+      registry.get<TH1>("hImpactParameterProngSqSum")->Fill(candidate.impactParameterProngSqSum());
+      registry.get<TH1>("hDecayLengthError")->Fill(candidate.errorDecayLength());
+      registry.get<TH1>("hDecayLengthXYError")->Fill(candidate.errorDecayLengthXY());
+      registry.get<TH1>("hImpactParameterError")->Fill(candidate.errorImpactParameter0());
+      registry.get<TH1>("hImpactParameterError")->Fill(candidate.errorImpactParameter1());
+      registry.get<TH1>("hImpactParameterError")->Fill(candidate.errorImpactParameter2());
+      registry.get<TH1>("hPtProng0")->Fill(candidate.ptProng0());
+      registry.get<TH1>("hPtProng1")->Fill(candidate.ptProng1());
+      registry.get<TH1>("hPtProng2")->Fill(candidate.ptProng2());
+      registry.get<TH1>("hd0Prong0")->Fill(candidate.impactParameter0());
+      registry.get<TH1>("hd0Prong1")->Fill(candidate.impactParameter1());
+      registry.get<TH1>("hd0Prong2")->Fill(candidate.impactParameter2());
+      //registry.get<TH1>("hSelectionStatus")->Fill(candidate.isSelDPlus());
     }
   }
 };