[PWGLF] Update decay3body analysis (#12317)

Author: creetz16

PWGLF/DataModel/Vtx3BodyTables.h

@@ -67,10 +67,10 @@ DECLARE_SOA_COLUMN(DCAZTrackPiToPV, dcaZtrackPiToPv, float);   //! DCAZ of pion
 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
+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(DCAVtxToDaughtersAv, dcaVtxToDaughtersAv, float); //! Quadratic sum of DCA between daughters at SV
 
 // CosPA
 DECLARE_SOA_COLUMN(CosPA, cosPA, float); //! Cosine of pointing angle of the 3body candidate
@@ -211,7 +211,7 @@ DECLARE_SOA_TABLE(Vtx3BodyDatas, "AOD", "VTX3BODYDATA", //!
                   vtx3body::DCAXYTrackPrToPV, vtx3body::DCAXYTrackPiToPV, vtx3body::DCAXYTrackDeToPV,
                   vtx3body::DCAZTrackPrToPV, vtx3body::DCAZTrackPiToPV, vtx3body::DCAZTrackDeToPV,
                   vtx3body::DCATrackPrToSV, vtx3body::DCATrackPiToSV, vtx3body::DCATrackDeToSV,
-                  vtx3body::DCAVtxDaughters,
+                  vtx3body::DCAVtxToDaughtersAv,
                   vtx3body::CosPA, vtx3body::Ct,
                   vtx3body::TPCNSigmaPr, vtx3body::TPCNSigmaPi, vtx3body::TPCNSigmaDe, vtx3body::TPCNSigmaPiBach,
                   vtx3body::TOFNSigmaDe,
@@ -260,7 +260,7 @@ DECLARE_SOA_TABLE(McVtx3BodyDatas, "AOD", "MC3BODYDATA", //!
                   vtx3body::DCAXYTrackPrToPV, vtx3body::DCAXYTrackPiToPV, vtx3body::DCAXYTrackDeToPV,
                   vtx3body::DCAZTrackPrToPV, vtx3body::DCAZTrackPiToPV, vtx3body::DCAZTrackDeToPV,
                   vtx3body::DCATrackPrToSV, vtx3body::DCATrackPiToSV, vtx3body::DCATrackDeToSV,
-                  vtx3body::DCAVtxDaughters,
+                  vtx3body::DCAVtxToDaughtersAv,
                   vtx3body::CosPA, vtx3body::Ct,
                   vtx3body::TPCNSigmaPr, vtx3body::TPCNSigmaPi, vtx3body::TPCNSigmaDe, vtx3body::TPCNSigmaPiBach,
                   vtx3body::TOFNSigmaDe,

PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx

@@ -143,7 +143,9 @@ struct decay3bodyBuilder {
     // building options
     Configurable<bool> useKFParticle{"useKFParticle", false, "Use KFParticle for decay3body building"};
     Configurable<bool> kfSetTopologicalConstraint{"kfSetTopologicalConstraint", false, "Set topological vertex constraint in case of KFParticle reconstruction"};
+    Configurable<bool> buildOnlyTracked{"buildOnlyTracked", false, "Build only tracked decay3bodys"};
     Configurable<bool> useSelections{"useSelections", true, "Apply selections during decay3body building"};
+    Configurable<bool> useChi2Selection{"useChi2Selection", true, "Apply chi2 selection during decay3body building"};
     Configurable<bool> useTPCforPion{"useTPCforPion", false, "Flag to ask for TPC info for pion track (PID, nClusters), false: pion track can be ITS only"};
     Configurable<bool> acceptTPCOnly{"acceptTPCOnly", false, "Accept TPC only tracks as daughters"};
     Configurable<bool> askOnlyITSMatch{"askOnlyITSMatch", true, "ask only ITS match to distinguish TPC only tracks"};
@@ -166,7 +168,7 @@ struct decay3bodyBuilder {
     Configurable<float> minTOFnSigmaDeuteron{"minTOFnSigmaDeuteron", -5.0, "Min TOF nSigma of deuteron daughter"};
     Configurable<float> maxTOFnSigmaDeuteron{"maxTOFnSigmaDeuteron", 5.0, "Max TOF nSigma of deuteron daughter"};
     Configurable<float> minPDeuteronUseTOF{"minPDeuteronUseTOF", 1.0, "Min P of deuteron to use TOF PID"};
-    Configurable<float> maxDCADauAtSV{"maxDCADauAtSV", 0.5, "Max DCA of daughters at SV (quadratic sum of daughter DCAs between each other)"};
+    Configurable<float> maxDCADauToSVaverage{"maxDCADauToSVaverage", 0.5, "Max DCA of daughters to SV (quadratic sum of daughter DCAs to SV / 3)"};
     // candidate selections
     Configurable<float> maxRapidity{"maxRapidity", 1.0, "Max rapidity of decay3body vertex"};
     Configurable<float> minPt{"minPt", 2.0, "Min Pt of decay3body candidate"};
@@ -326,7 +328,7 @@ struct decay3bodyBuilder {
     helper.decay3bodyselections.minTOFnSigmaDeuteron = decay3bodyBuilderOpts.minTOFnSigmaDeuteron;
     helper.decay3bodyselections.maxTOFnSigmaDeuteron = decay3bodyBuilderOpts.maxTOFnSigmaDeuteron;
     helper.decay3bodyselections.minPDeuteronUseTOF = decay3bodyBuilderOpts.minPDeuteronUseTOF;
-    helper.decay3bodyselections.maxDCADauAtSV = decay3bodyBuilderOpts.maxDCADauAtSV;
+    helper.decay3bodyselections.maxDCADauToSVaverage = decay3bodyBuilderOpts.maxDCADauToSVaverage;
     helper.decay3bodyselections.maxRapidity = decay3bodyBuilderOpts.maxRapidity;
     helper.decay3bodyselections.minPt = decay3bodyBuilderOpts.minPt;
     helper.decay3bodyselections.maxPt = decay3bodyBuilderOpts.maxPt;
@@ -393,7 +395,7 @@ struct decay3bodyBuilder {
     LOGF(info, "-~> min TOF nSigma deuteron ......: %f", decay3bodyBuilderOpts.minTOFnSigmaDeuteron.value);
     LOGF(info, "-~> max TOF nSigma deuteron ......: %f", decay3bodyBuilderOpts.maxTOFnSigmaDeuteron.value);
     LOGF(info, "-~> min p bach use TOF ...........: %f", decay3bodyBuilderOpts.minPDeuteronUseTOF.value);
-    LOGF(info, "-~> max DCA dau at SV ............: %f", decay3bodyBuilderOpts.maxDCADauAtSV.value);
+    LOGF(info, "-~> max DCA dau at SV ............: %f", decay3bodyBuilderOpts.maxDCADauToSVaverage.value);
     LOGF(info, "-~> max rapidity .................: %f", decay3bodyBuilderOpts.maxRapidity.value);
     LOGF(info, "-~> min pT .......................: %f", decay3bodyBuilderOpts.minPt.value);
     LOGF(info, "-~> max pT .......................: %f", decay3bodyBuilderOpts.maxPt.value);
@@ -720,6 +722,11 @@ struct decay3bodyBuilder {
     registry.fill(HIST("Counters/hInputStatistics"), kVtx3BodyDatas, decay3bodys.size());
     int lastRunNumber = -1;
     for (const auto& decay3body : decay3bodys) {
+      // only build tracked decay3body if aksed
+      if (decay3bodyBuilderOpts.buildOnlyTracked && fTrackedClSizeVector[decay3body.globalIndex()] == 0) {
+        continue;
+      }
+
       // skip decay3body without assigned collision
       /// TODO: do we want this??
       if (decay3body.collisionId() < 0) {
@@ -782,6 +789,7 @@ struct decay3bodyBuilder {
                                            decay3bodyBuilderOpts.useKFParticle,
                                            decay3bodyBuilderOpts.kfSetTopologicalConstraint,
                                            decay3bodyBuilderOpts.useSelections,
+                                           decay3bodyBuilderOpts.useChi2Selection,
                                            decay3bodyBuilderOpts.useTPCforPion,
                                            decay3bodyBuilderOpts.acceptTPCOnly,
                                            decay3bodyBuilderOpts.askOnlyITSMatch,
@@ -958,7 +966,7 @@ struct decay3bodyBuilder {
                                    -1., -1., -1.,      // trackDCAxyToPV: 0 - proton, 1 - pion, 2 - deuteron
                                    -1., -1., -1.,      // trackDCAzToPV: 0 - proton, 1 - pion, 2 - deuteron
                                    -1., -1., -1.,      // daughterDCAtoSV: 0 - proton, 1 - pion, 2 - deuteron
-                                   -1.,                // daughterDCAatSV
+                                   -1.,                // daughterDCAtoSVaverage
                                    -1., -1.,           // cosPA, ctau
                                    -1., -1., -1., -1., // tpcNsigma: 0 - proton, 1 - pion, 2 - deuteron, 3 - bach with pion hyp
                                    -1.,                // tofNsigmaDeuteron
@@ -1117,7 +1125,7 @@ struct decay3bodyBuilder {
                              helper.decay3body.trackDCAxyToPV[0], helper.decay3body.trackDCAxyToPV[1], helper.decay3body.trackDCAxyToPV[2],    // 0 - proton, 1 - pion, 2 - deuteron
                              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.daughterDCAtoSVaverage,
                              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,
@@ -1146,7 +1154,7 @@ struct decay3bodyBuilder {
                                helper.decay3body.trackDCAxyToPV[0], helper.decay3body.trackDCAxyToPV[1], helper.decay3body.trackDCAxyToPV[2],    // 0 - proton, 1 - pion, 2 - deuteron
                                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.daughterDCAtoSVaverage,
                                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,
@@ -1184,6 +1192,7 @@ struct decay3bodyBuilder {
                                         decay3bodyBuilderOpts.useKFParticle,
                                         decay3bodyBuilderOpts.kfSetTopologicalConstraint,
                                         decay3bodyBuilderOpts.useSelections,
+                                        decay3bodyBuilderOpts.useChi2Selection,
                                         decay3bodyBuilderOpts.useTPCforPion,
                                         decay3bodyBuilderOpts.acceptTPCOnly,
                                         decay3bodyBuilderOpts.askOnlyITSMatch,

PWGLF/Utils/decay3bodyBuilderHelper.h

@@ -74,15 +74,11 @@ struct decay3bodyCandidate {
   int sign;
   float momentum[3];
   float position[3];
-  // std::array<float, 3> momentum = {0.0f, 0.0f, 0.0f};
-  // std::array<float, 3> position = {0.0f, 0.0f, 0.0f};
-  // float dcaToPV = 0.0f;
-  // 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
+  float daughterDCAtoSVaverage = 0.0f;                       // average of quadratic sum of daughter DCAs to SV
   std::array<float, 3> daughterDCAtoSV = {0.0f, 0.0f, 0.0f}; // 0 - pos, 1 - neg, 2 - bach
 
   // covariance matrix
@@ -148,7 +144,7 @@ class decay3bodyBuilderHelper
     double minTOFnSigmaDeuteron;
     double maxTOFnSigmaDeuteron;
     float minPDeuteronUseTOF;
-    float maxDCADauAtSV;
+    float maxDCADauToSVaverage;
     // candidate
     float maxRapidity;
     float minPt;
@@ -188,6 +184,7 @@ class decay3bodyBuilderHelper
                                 bool useKFParticle = false,
                                 bool kfSetTopologicalConstraint = false,
                                 bool useSelections = true,
+                                bool useChi2Selection = true,
                                 bool useTPCforPion = false,
                                 bool acceptTPCOnly = false,
                                 bool askOnlyITSMatch = true,
@@ -342,18 +339,6 @@ class decay3bodyBuilderHelper
     auto trackPionPt = std::sqrt(decay3body.momPion[0] * decay3body.momPion[0] + decay3body.momPion[1] * decay3body.momPion[1]);
     auto trackDeuteronPt = std::sqrt(decay3body.momDeuteron[0] * decay3body.momDeuteron[0] + decay3body.momDeuteron[1] * decay3body.momDeuteron[1]);
 
-    // DCA between daughters at SV
-    decay3body.daughterDCAatSV = std::hypot(
-      std::hypot(decay3body.posProton[0] - decay3body.posPion[0],
-                 decay3body.posProton[1] - decay3body.posPion[1],
-                 decay3body.posProton[2] - decay3body.posPion[2]),
-      std::hypot(decay3body.posProton[0] - decay3body.posDeuteron[0],
-                 decay3body.posProton[1] - decay3body.posDeuteron[1],
-                 decay3body.posProton[2] - decay3body.posDeuteron[2]),
-      std::hypot(decay3body.posPion[0] - decay3body.posProton[0],
-                 decay3body.posPion[1] - decay3body.posProton[1],
-                 decay3body.posPion[2] - decay3body.posProton[2]));
-
     // daughter DCA to SV
     // proton daughter
     decay3body.daughterDCAtoSV[0] = std::hypot(
@@ -371,6 +356,12 @@ class decay3bodyBuilderHelper
       decay3body.posDeuteron[1] - decay3body.position[1],
       decay3body.posDeuteron[2] - decay3body.position[2]);
 
+    // DCA daughters to SV average of quadratic sum
+    decay3body.daughterDCAtoSVaverage = (decay3body.daughterDCAtoSV[0] * decay3body.daughterDCAtoSV[0] +
+                                         decay3body.daughterDCAtoSV[1] * decay3body.daughterDCAtoSV[1] +
+                                         decay3body.daughterDCAtoSV[2] * decay3body.daughterDCAtoSV[2]) /
+                                        3;
+
     //_____________________________________________________
     // selections after vertex fit
     if (useSelections) {
@@ -392,7 +383,7 @@ class decay3bodyBuilderHelper
       }
 
       // daughter DCAs at SV
-      if (decay3body.daughterDCAatSV > decay3bodyselections.maxDCADauAtSV) {
+      if (decay3body.daughterDCAtoSVaverage > decay3bodyselections.maxDCADauToSVaverage) {
         decay3body = {};
         return false;
       }
@@ -418,7 +409,7 @@ class decay3bodyBuilderHelper
       }
 
       // vertex chi2
-      if (decay3body.chi2 > decay3bodyselections.maxChi2) {
+      if (useChi2Selection && decay3body.chi2 > decay3bodyselections.maxChi2) {
         decay3body = {};
         return false;
       }