[PWGLF] Store cosPA and ct in datamodel for h3l_3body

PWGLF/DataModel/Vtx3BodyTables.h

@@ -17,10 +17,12 @@
 #ifndef PWGLF_DATAMODEL_VTX3BODYTABLES_H_
 #define PWGLF_DATAMODEL_VTX3BODYTABLES_H_
 
-#include <cmath>
-#include "Framework/AnalysisDataModel.h"
 #include "Common/Core/RecoDecay.h"
+
 #include "CommonConstants/PhysicsConstants.h"
+#include "Framework/AnalysisDataModel.h"
+
+#include <cmath>
 
 namespace o2::aod
 {
@@ -30,8 +32,8 @@
 DECLARE_SOA_INDEX_COLUMN_FULL(TrackPr, trackPr, int, Tracks, "_pr"); //!
 DECLARE_SOA_INDEX_COLUMN_FULL(TrackPi, trackPi, int, Tracks, "_pi"); //!
 DECLARE_SOA_INDEX_COLUMN_FULL(TrackDe, trackDe, int, Tracks, "_de"); //!
-DECLARE_SOA_INDEX_COLUMN(Collision, collision);                   //!
-DECLARE_SOA_INDEX_COLUMN(Decay3Body, decay3body);                 //!
+DECLARE_SOA_INDEX_COLUMN(Collision, collision);                      //!
+DECLARE_SOA_INDEX_COLUMN(Decay3Body, decay3body);                    //!
 
 // General 3 body Vtx properties
 DECLARE_SOA_COLUMN(Mass, mass, float); //! candidate mass (with H3L or Anti-H3L mass hypothesis depending on deuteron charge)
@@ -57,19 +59,25 @@
 DECLARE_SOA_COLUMN(PzTrackDe, pzTrackDe, float); //! track2 pz at min
 
 // DCAs to PV
 DECLARE_SOA_COLUMN(DCAXYTrackPrToPV, dcaXYtrackPrToPv, float); //! DCAXY of proton to PV
 DECLARE_SOA_COLUMN(DCAXYTrackPiToPV, dcaXYtrackPiToPv, float); //! DCAXY of pion to PV
 DECLARE_SOA_COLUMN(DCAXYTrackDeToPV, dcaXYtrackDeToPv, float); //! DCAXY of deuteron to PV
 DECLARE_SOA_COLUMN(DCAZTrackPrToPV, dcaZtrackPrToPv, float);   //! DCAZ of proton to PV
 DECLARE_SOA_COLUMN(DCAZTrackPiToPV, dcaZtrackPiToPv, float);   //! DCAZ of pion to PV
 DECLARE_SOA_COLUMN(DCAZTrackDeToPV, dcaZtrackDeToPv, float);   //! DCAZ of deuteron to PV
 
 // DCAs to SV
 DECLARE_SOA_COLUMN(DCATrackPrToSV, dcaTrackPrToSv, float);   //! DCA of proton to SV
 DECLARE_SOA_COLUMN(DCATrackPiToSV, dcaTrackPiToSv, float);   //! DCA of pion to SV
 DECLARE_SOA_COLUMN(DCATrackDeToSV, dcaTrackDeToSv, float);   //! DCA of deuteron to SV
 DECLARE_SOA_COLUMN(DCAVtxDaughters, dcaVtxdaughters, float); //! Quadratic sum of DCA between daughters at SV
 
+// CosPA
+DECLARE_SOA_COLUMN(CosPA, cosPA, float); //! Cosine of pointing angle of the 3body candidate
+
+// Ct
+DECLARE_SOA_COLUMN(Ct, ct, float); //! Reconstruction Ct of 3body candidate
+
 // Strangeness tracking
 DECLARE_SOA_COLUMN(TrackedClSize, trackedClSize, float); //! Average ITS cluster size of strangeness tracked 3body
 
@@ -141,14 +149,10 @@
 // Distance Over To Mom
 DECLARE_SOA_DYNAMIC_COLUMN(DistOverTotMom, distovertotmom, //! PV to 3 body decay distance over total momentum
                            [](float X, float Y, float Z, float Px, float Py, float Pz, float pvX, float pvY, float pvZ) {
                              float P = RecoDecay::sqrtSumOfSquares(Px, Py, Pz);
                              return std::sqrt(std::pow(X - pvX, 2) + std::pow(Y - pvY, 2) + std::pow(Z - pvZ, 2)) / (P + 1E-10);
                            });
 
-// CosPA
-DECLARE_SOA_DYNAMIC_COLUMN(VtxCosPA, vtxcosPA, //! 3 body vtx CosPA
-                           [](float X, float Y, float Z, float Px, float Py, float Pz, float pvX, float pvY, float pvZ) -> float { return RecoDecay::cpa(std::array{pvX, pvY, pvZ}, std::array{X, Y, Z}, std::array{Px, Py, Pz}); });
-
 // Dca to PV
 DECLARE_SOA_DYNAMIC_COLUMN(DCAVtxToPV, dcavtxtopv, //! DCA of 3 body vtx to PV
                            [](float X, float Y, float Z, float Px, float Py, float Pz, float pvX, float pvY, float pvZ) -> float { return std::sqrt((std::pow((pvY - Y) * Pz - (pvZ - Z) * Py, 2) + std::pow((pvX - X) * Pz - (pvZ - Z) * Px, 2) + std::pow((pvX - X) * Py - (pvY - Y) * Px, 2)) / (Px * Px + Py * Py + Pz * Pz)); });
@@ -207,6 +211,7 @@
                   vtx3body::DCAZTrackPrToPV, vtx3body::DCAZTrackPiToPV, vtx3body::DCAZTrackDeToPV,
                   vtx3body::DCATrackPrToSV, vtx3body::DCATrackPiToSV, vtx3body::DCATrackDeToSV,
                   vtx3body::DCAVtxDaughters,
+                  vtx3body::CosPA, vtx3body::Ct,
                   vtx3body::TPCNSigmaPr, vtx3body::TPCNSigmaPi, vtx3body::TPCNSigmaDe, vtx3body::TPCNSigmaPiBach,
                   vtx3body::TOFNSigmaDe,
                   vtx3body::ITSClSizePr, vtx3body::ITSClSizePi, vtx3body::ITSClSizeDe,
@@ -218,7 +223,6 @@
                   vtx3body::Pt<vtx3body::Px, vtx3body::Py>,
                   vtx3body::VtxRadius<vtx3body::X, vtx3body::Y>,
                   vtx3body::DistOverTotMom<vtx3body::X, vtx3body::Y, vtx3body::Z, vtx3body::Px, vtx3body::Py, vtx3body::Pz>,
-                  vtx3body::VtxCosPA<vtx3body::X, vtx3body::Y, vtx3body::Z, vtx3body::Px, vtx3body::Py, vtx3body::Pz>,
                   vtx3body::DCAVtxToPV<vtx3body::X, vtx3body::Y, vtx3body::Z, vtx3body::Px, vtx3body::Py, vtx3body::Pz>,
 
                   // Longitudinal
@@ -256,6 +260,7 @@
                   vtx3body::DCAZTrackPrToPV, vtx3body::DCAZTrackPiToPV, vtx3body::DCAZTrackDeToPV,
                   vtx3body::DCATrackPrToSV, vtx3body::DCATrackPiToSV, vtx3body::DCATrackDeToSV,
                   vtx3body::DCAVtxDaughters,
+                  vtx3body::CosPA, vtx3body::Ct,
                   vtx3body::TPCNSigmaPr, vtx3body::TPCNSigmaPi, vtx3body::TPCNSigmaDe, vtx3body::TPCNSigmaPiBach,
                   vtx3body::TOFNSigmaDe,
                   vtx3body::ITSClSizePr, vtx3body::ITSClSizePi, vtx3body::ITSClSizeDe,
@@ -283,7 +288,6 @@
                   vtx3body::VtxRadius<vtx3body::X, vtx3body::Y>,
                   vtx3body::GenRadius<vtx3body::GenX, vtx3body::GenY>,
                   vtx3body::DistOverTotMom<vtx3body::X, vtx3body::Y, vtx3body::Z, vtx3body::Px, vtx3body::Py, vtx3body::Pz>,
-                  vtx3body::VtxCosPA<vtx3body::X, vtx3body::Y, vtx3body::Z, vtx3body::Px, vtx3body::Py, vtx3body::Pz>,
                   vtx3body::DCAVtxToPV<vtx3body::X, vtx3body::Y, vtx3body::Z, vtx3body::Px, vtx3body::Py, vtx3body::Pz>,
 
                   // Longitudinal

PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx

@@ -958,6 +958,7 @@ struct decay3bodyBuilder {
                                    -1., -1., -1.,      // trackDCAzToPV: 0 - proton, 1 - pion, 2 - deuteron
                                    -1., -1., -1.,      // daughterDCAtoSV: 0 - proton, 1 - pion, 2 - deuteron
                                    -1.,                // daughterDCAatSV
+                                   -1., -1.,           // cosPA, ctau
                                    -1., -1., -1., -1., // tpcNsigma: 0 - proton, 1 - pion, 2 - deuteron, 3 - bach with pion hyp
                                    -1.,                // tofNsigmaDeuteron
                                    -1., -1., -1.,      // average ITS cluster sizes: proton, pion, deuteron
@@ -1115,6 +1116,7 @@ struct decay3bodyBuilder {
                              helper.decay3body.trackDCAzToPV[0], helper.decay3body.trackDCAzToPV[1], helper.decay3body.trackDCAzToPV[2],       // 0 - proton, 1 - pion, 2 - deuteron
                              helper.decay3body.daughterDCAtoSV[0], helper.decay3body.daughterDCAtoSV[1], helper.decay3body.daughterDCAtoSV[2], // 0 - proton, 1 - pion, 2 - deuteron
                              helper.decay3body.daughterDCAatSV,
+                             helper.decay3body.cosPA, helper.decay3body.ctau,
                              helper.decay3body.tpcNsigma[0], helper.decay3body.tpcNsigma[1], helper.decay3body.tpcNsigma[2], helper.decay3body.tpcNsigma[2], // 0 - proton, 1 - pion, 2 - deuteron, 3 - bach with pion hyp
                              helper.decay3body.tofNsigmaDeuteron,
                              helper.decay3body.averageITSClSize[0], helper.decay3body.averageITSClSize[1], helper.decay3body.averageITSClSize[2], // 0 - proton, 1 - pion, 2 - deuteron
@@ -1143,6 +1145,7 @@ struct decay3bodyBuilder {
                                helper.decay3body.trackDCAzToPV[0], helper.decay3body.trackDCAzToPV[1], helper.decay3body.trackDCAzToPV[2],       // 0 - proton, 1 - pion, 2 - deuteron
                                helper.decay3body.daughterDCAtoSV[0], helper.decay3body.daughterDCAtoSV[1], helper.decay3body.daughterDCAtoSV[2], // 0 - proton, 1 - pion, 2 - deuteron
                                helper.decay3body.daughterDCAatSV,
+                               helper.decay3body.cosPA, helper.decay3body.ctau,
                                helper.decay3body.tpcNsigma[0], helper.decay3body.tpcNsigma[1], helper.decay3body.tpcNsigma[2], helper.decay3body.tpcNsigma[2], // 0 - proton, 1 - pion, 2 - deuteron, 3 - bach with pion hyp
                                helper.decay3body.tofNsigmaDeuteron,
                                helper.decay3body.averageITSClSize[0], helper.decay3body.averageITSClSize[1], helper.decay3body.averageITSClSize[2], // 0 - proton, 1 - pion, 2 - deuteron

PWGLF/Utils/decay3bodyBuilderHelper.h

@@ -12,27 +12,29 @@
 #ifndef PWGLF_UTILS_DECAY3BODYBUILDERHELPER_H_
 #define PWGLF_UTILS_DECAY3BODYBUILDERHELPER_H_
 
-#include <cstdlib>
-#include <cmath>
-#include <array>
+#include "Common/Core/RecoDecay.h"
+#include "Common/Core/trackUtilities.h"
+#include "Tools/KFparticle/KFUtilities.h"
+
+#include "CommonConstants/PhysicsConstants.h"
 #include "DCAFitter/DCAFitterN.h"
-#include "Framework/AnalysisDataModel.h"
-#include "ReconstructionDataFormats/Track.h"
 #include "DetectorsBase/GeometryManager.h"
 #include "DetectorsVertexing/SVertexHypothesis.h"
-#include "CommonConstants/PhysicsConstants.h"
-#include "Common/Core/trackUtilities.h"
-#include "Common/Core/RecoDecay.h"
-#include "Tools/KFparticle/KFUtilities.h"
+#include "Framework/AnalysisDataModel.h"
+#include "ReconstructionDataFormats/Track.h"
+
+#include <array>
+#include <cmath>
+#include <cstdlib>
 
 #ifndef HomogeneousField
 #define HomogeneousField
 #endif
 
 /// includes KFParticle
-#include "KFParticle.h"
 #include "KFPTrack.h"
 #include "KFPVertex.h"
+#include "KFParticle.h"
 #include "KFParticleBase.h"
 #include "KFVertex.h"
 
@@ -78,6 +80,8 @@ struct decay3bodyCandidate {
   // float dcaxyToPV = 0.0f;
   float chi2 = 0.0f;
   float trackedClSize = 0.0f;
+  float cosPA = 0.0f;                                        // cosine of pointing angle
+  float ctau = 0.0f;                                         // ctau of the candidate
   float daughterDCAatSV = 0.0f;                              // quadratic sum of DCA between daughters at SV
   std::array<float, 3> daughterDCAtoSV = {0.0f, 0.0f, 0.0f}; // 0 - pos, 1 - neg, 2 - bach
 
@@ -413,27 +417,33 @@ class decay3bodyBuilderHelper
         return false;
       }
 
-      // pointing angle
-      float cpa = RecoDecay::cpa(std::array{pvX, pvY, pvZ}, std::array{decay3body.position[0], decay3body.position[1], decay3body.position[2]}, std::array{decay3body.momentum[0], decay3body.momentum[1], decay3body.momentum[2]});
-      if (cpa < decay3bodyselections.minCosPA) {
+      // vertex chi2
+      if (decay3body.chi2 > decay3bodyselections.maxChi2) {
         decay3body = {};
         return false;
       }
+    }
 
-      // vertex chi2
-      if (decay3body.chi2 > decay3bodyselections.maxChi2) {
+    // pointing angle
+    float cpa = RecoDecay::cpa(std::array{pvX, pvY, pvZ}, std::array{decay3body.position[0], decay3body.position[1], decay3body.position[2]}, std::array{decay3body.momentum[0], decay3body.momentum[1], decay3body.momentum[2]});
+    if (useSelections) {
+      if (cpa < decay3bodyselections.minCosPA) {
         decay3body = {};
         return false;
       }
+    }
+    decay3body.cosPA = cpa;
 
-      // ctau
-      float P = RecoDecay::sqrtSumOfSquares(decay3body.momentum[0], decay3body.momentum[1], decay3body.momentum[2]);
-      float ctau = std::sqrt(std::pow(decay3body.position[0] - pvX, 2) + std::pow(decay3body.position[1] - pvY, 2) + std::pow(decay3body.position[2] - pvZ, 2)) / (P + 1E-10) * o2::constants::physics::MassHyperTriton;
+    // ctau
+    float P = RecoDecay::sqrtSumOfSquares(decay3body.momentum[0], decay3body.momentum[1], decay3body.momentum[2]);
+    float ctau = std::sqrt(std::pow(decay3body.position[0] - pvX, 2) + std::pow(decay3body.position[1] - pvY, 2) + std::pow(decay3body.position[2] - pvZ, 2)) / (P + 1E-10) * o2::constants::physics::MassHyperTriton;
+    if (useSelections) {
       if (ctau < decay3bodyselections.minCtau || ctau > decay3bodyselections.maxCtau) {
         decay3body = {};
         return false;
       }
     }
+    decay3body.ctau = ctau;
 
     //_______________________________________________________________________
     // SVertexer selections in case of event mixing