diff --git a/PWGLF/DataModel/ReducedHeptaQuarkTables.h b/PWGLF/DataModel/ReducedHeptaQuarkTables.h
index b547933cb59..76cc71beb09 100644
--- a/PWGLF/DataModel/ReducedHeptaQuarkTables.h
+++ b/PWGLF/DataModel/ReducedHeptaQuarkTables.h
@@ -38,6 +38,8 @@ DECLARE_SOA_TABLE(RedHQEvents, "AOD", "REDHQEVENT",
                   bc::GlobalBC,
                   bc::RunNumber,
                   timestamp::Timestamp,
+                  collision::PosX,
+                  collision::PosY,
                   collision::PosZ,
                   collision::NumContrib,
                   redhqevent::Centrality,
diff --git a/PWGLF/TableProducer/Resonances/HeptaQuarktable.cxx b/PWGLF/TableProducer/Resonances/HeptaQuarktable.cxx
index 818c43da234..1d6b20dd6e4 100644
--- a/PWGLF/TableProducer/Resonances/HeptaQuarktable.cxx
+++ b/PWGLF/TableProducer/Resonances/HeptaQuarktable.cxx
@@ -438,7 +438,7 @@ struct heptaquarktable {
     if (keepEventDoubleHQ && numberPhi > 1 && numberLambda > 0 && (hqresonance.size() == hqresonanced1.size()) && (hqresonance.size() == hqresonanced2.size())) {
       histos.fill(HIST("hEventstat"), 2.5);
       /////////// Fill collision table///////////////
-      redHQEvents(bc.globalBC(), currentRunNumber, bc.timestamp(), collision.posZ(), collision.numContrib(), centrality, numberPhi, numberLambda);
+      redHQEvents(bc.globalBC(), currentRunNumber, bc.timestamp(), collision.posX(), collision.posY(), collision.posZ(), collision.numContrib(), centrality, numberPhi, numberLambda);
       auto indexEvent = redHQEvents.lastIndex();
       //// Fill track table for HQ//////////////////
       for (auto if1 = hqresonance.begin(); if1 != hqresonance.end(); ++if1) {
diff --git a/PWGLF/Tasks/Resonances/heptaquark.cxx b/PWGLF/Tasks/Resonances/heptaquark.cxx
index 407c3199073..f7c30a46ddf 100644
--- a/PWGLF/Tasks/Resonances/heptaquark.cxx
+++ b/PWGLF/Tasks/Resonances/heptaquark.cxx
@@ -12,28 +12,33 @@
 /// \brief this is a starting point for the Resonances tutorial
 /// \author junlee kim
 
+#include "PWGLF/DataModel/ReducedHeptaQuarkTables.h"
+
+#include "Common/Core/trackUtilities.h"
+
+#include "CommonConstants/PhysicsConstants.h"
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/StepTHn.h"
+#include "Framework/runDataProcessing.h"
 #include <Framework/Configurable.h>
-#include <TLorentzVector.h>
+
 #include <Math/GenVector/Boost.h>
-#include <Math/Vector4D.h>
 #include <Math/Vector3D.h>
+#include <Math/Vector4D.h>
+#include <TLorentzVector.h>
 #include <TMath.h>
 #include <TRandom3.h>
+
 #include <fairlogger/Logger.h>
+
+#include <algorithm>
 #include <iostream>
 #include <iterator>
 #include <string>
 #include <vector>
 
-#include "Framework/AnalysisTask.h"
-#include "Framework/ASoAHelpers.h"
-#include "Framework/runDataProcessing.h"
-#include "Framework/AnalysisDataModel.h"
-#include "Framework/StepTHn.h"
-#include "Common/Core/trackUtilities.h"
-#include "PWGLF/DataModel/ReducedHeptaQuarkTables.h"
-#include "CommonConstants/PhysicsConstants.h"
-
 using namespace o2;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
@@ -47,14 +52,22 @@ struct heptaquark {
   Configurable<int> cfgPIDStrategy{"cfgPIDStrategy", 3, "PID strategy 1"};
   Configurable<float> cfgPIDPrPi{"cfgPIDPrPi", 3, "PID selection for proton and pion"};
 
-  Configurable<float> minPhiMass{"minPhiMass", 1.01, "Minimum phi mass"};
-  Configurable<float> maxPhiMass{"maxPhiMass", 1.03, "Maximum phi mass"};
+  Configurable<float> cfgMinPhiMass{"cfgMinPhiMass", 1.01, "Minimum phi mass"};
+  Configurable<float> cfgMaxPhiMass{"cfgMaxPhiMass", 1.03, "Maximum phi mass"};
 
-  Configurable<float> minLambdaMass{"minLambdaMass", 1.1, "Minimum lambda mass"};
-  Configurable<float> maxLambdaMass{"maxLambdaMass", 1.13, "Maximum lambda mass"};
+  Configurable<float> cfgMinLambdaMass{"cfgMinLambdaMass", 1.1, "Minimum lambda mass"};
+  Configurable<float> cfgMaxLambdaMass{"cfgMaxLambdaMass", 1.13, "Maximum lambda mass"};
 
-  Configurable<float> cutNsigmaTPC{"cutNsigmaTPC", 2.5, "nsigma cut TPC"};
-  Configurable<float> cutNsigmaTOF{"cutNsigmaTOF", 3.0, "nsigma cut TOF"};
+  Configurable<float> cfgNsigmaTPC{"cfgNsigmaTPC", 2.5, "nsigma cut TPC"};
+  Configurable<float> cfgNsigmaTOF{"cfgNsigmaTOF", 3.0, "nsigma cut TOF"};
+
+  Configurable<bool> cfgSelectHQ{"cfgSelectHQ", true, "switch to select HQ"};
+
+  Configurable<float> cfgMinPhiPt{"cfgMinPhiPt", 0.2, "Minimum phi pt"};
+  Configurable<float> cfgMinLambdaPt{"cfgMinLambdaPt", 0.5, "Minimum lambda pt"};
+
+  Configurable<float> cfgSoftFraction{"cfgSoftFraction", 0.01, "Minimum allowed softest fraction"};
+  Configurable<float> cfgCollinear{"cfgCollinear", 0.98, "Maximum allowed collinear selection"};
 
   ConfigurableAxis massAxis{"massAxis", {600, 2.8, 3.4}, "Invariant mass axis"};
   ConfigurableAxis ptAxis{"ptAxis", {VARIABLE_WIDTH, 0.2, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.5, 8.0, 10.0, 100.0}, "Transverse momentum bins"};
@@ -95,76 +108,76 @@ struct heptaquark {
   {
     if (PIDStrategy == 0) {
       if (TOFHit != 1) {
-        if (TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
+        if (TMath::Abs(nsigmaTPC) < cfgNsigmaTPC) {
           return true;
         }
       }
       if (TOFHit == 1) {
-        if (TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
+        if (TMath::Abs(nsigmaTOF) < cfgNsigmaTOF) {
           return true;
         }
       }
     }
     if (PIDStrategy == 1) {
       if (ptcand < 0.5) {
-        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
+        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cfgNsigmaTPC) {
           return true;
         }
-        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
+        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cfgNsigmaTOF) {
           return true;
         }
       }
       if (ptcand >= 0.5) {
-        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
+        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cfgNsigmaTOF) {
           return true;
         }
       }
     }
     if (PIDStrategy == 2) {
       if (ptcand < 0.5) {
-        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
+        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cfgNsigmaTPC) {
           return true;
         }
-        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
+        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cfgNsigmaTOF) {
           return true;
         }
       }
       if (ptcand >= 0.5 && ptcand < 1.2) {
-        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
+        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cfgNsigmaTOF) {
           return true;
         }
-        if (TOFHit != 1 && nsigmaTPC > -1.5 && nsigmaTPC < cutNsigmaTPC) {
+        if (TOFHit != 1 && nsigmaTPC > -1.5 && nsigmaTPC < cfgNsigmaTPC) {
           return true;
         }
       }
       if (ptcand >= 1.2) {
-        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
+        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cfgNsigmaTOF) {
           return true;
         }
-        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
+        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cfgNsigmaTPC) {
           return true;
         }
       }
     }
     if (PIDStrategy == 3) {
       if (ptcand < 0.5) {
-        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
+        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cfgNsigmaTPC) {
           return true;
         }
-        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
+        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cfgNsigmaTOF) {
           return true;
         }
       }
       if (ptcand >= 0.5 && ptcand < 1.2) {
-        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
+        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cfgNsigmaTOF) {
           return true;
         }
       }
       if (ptcand >= 1.2) {
-        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
+        if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cfgNsigmaTOF) {
           return true;
         }
-        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
+        if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cfgNsigmaTPC) {
           return true;
         }
       }
@@ -172,28 +185,38 @@ struct heptaquark {
     return false;
   }
 
-  template <typename V01, typename V02>
-  ROOT::Math::XYZVector getDCAofV0V0(V01 const& v01, V02 const& v02)
-  {
-    ROOT::Math::XYZVector v01pos, v02pos, v01mom, v02mom;
-    v01pos.SetXYZ(v01.x(), v01.y(), v01.z());
-    v02pos.SetXYZ(v02.x(), v02.y(), v02.z());
-    v01mom.SetXYZ(v01.px(), v01.py(), v01.pz());
-    v02mom.SetXYZ(v02.px(), v02.py(), v02.pz());
-
-    ROOT::Math::XYZVector posdiff = v02pos - v01pos;
-    ROOT::Math::XYZVector cross = v01mom.Cross(v02mom);
-    ROOT::Math::XYZVector dcaVec = (posdiff.Dot(cross) / cross.Mag2()) * cross;
-    return dcaVec;
-  }
-
-  template <typename V01, typename V02>
-  float getCPA(V01 const& v01, V02 const& v02)
+  template <typename HQ1, typename HQ2, typename HQ3>
+  int selectHQ(HQ1 const& hq1, HQ2 const& hq2, HQ3 const& hq3)
   {
-    ROOT::Math::XYZVector v01mom, v02mom;
-    v01mom.SetXYZ(v01.px() / v01.p(), v01.py() / v01.p(), v01.pz() / v01.p());
-    v02mom.SetXYZ(v02.px() / v02.p(), v02.py() / v02.p(), v02.pz() / v02.p());
-    return v01mom.Dot(v02mom);
+    int selection = 0;
+    if (hq1.Pt() < cfgMinPhiPt || hq2.Pt() < cfgMinPhiPt || hq3.Pt() < cfgMinLambdaPt)
+      selection += 1;
+
+    double sumE = hq1.E() + hq2.E() + hq3.E();
+    double emin = std::min({hq1.E(), hq2.E(), hq3.E()});
+    double fmin = emin / std::max(1e-9, sumE);
+    if (fmin < cfgSoftFraction)
+      selection += 2;
+
+    auto ex = hq1 + hq2 + hq3;
+    TVector3 boost = -ex.BoostVector();
+    auto hqphipair_boost = hq1 + hq2;
+    auto hqlambda_boost = hq3;
+    hqphipair_boost.Boost(boost);
+    hqlambda_boost.Boost(boost);
+    double cosHel = hqlambda_boost.Vect().Dot(hqphipair_boost.Vect()) / (hqlambda_boost.Vect().Mag() * hqphipair_boost.Vect().Mag());
+    if (std::abs(cosHel) > cfgCollinear)
+      selection += 4;
+    /*
+        ROOT::Math::XYZVector rPV(col.posX(), col.posY(), col.posZ());
+        ROOT::Math::XYZVector rSV(hq3.hqx(), hq3.hqy(), hq3.hqz());
+        ROOT::Math::XYZVector L = rSV - rPV;
+        ROOT::Math::XYZVector exMom(ex.Px(), ex.Py(), ex.Pz());
+        double cosPoint = L.Dot(exMom) / (L.R() * pEx.R() + 1e-9);
+        if (cosPoint < cfgCosPoint)
+          return 8;
+    */
+    return selection;
   }
 
   ROOT::Math::PxPyPzMVector DauVec1, DauVec2;
@@ -214,7 +237,7 @@ struct heptaquark {
       if (hqtrackd1.hqId() != 333)
         continue;
 
-      if (hqtrackd1.hqMass() < minPhiMass || hqtrackd1.hqMass() > maxPhiMass)
+      if (hqtrackd1.hqMass() < cfgMinPhiMass || hqtrackd1.hqMass() > cfgMaxPhiMass)
         continue;
 
       DauVec1 = ROOT::Math::PxPyPzMVector(hqtrackd1.hqd1Px(), hqtrackd1.hqd1Py(), hqtrackd1.hqd1Pz(), massKa);
@@ -246,7 +269,7 @@ struct heptaquark {
         if (hqtrackd2.hqId() != 333)
           continue;
 
-        if (hqtrackd2.hqMass() < minPhiMass || hqtrackd2.hqMass() > maxPhiMass)
+        if (hqtrackd2.hqMass() < cfgMinPhiMass || hqtrackd2.hqMass() > cfgMaxPhiMass)
           continue;
 
         DauVec1 = ROOT::Math::PxPyPzMVector(hqtrackd2.hqd1Px(), hqtrackd2.hqd1Py(), hqtrackd2.hqd1Pz(), massKa);
@@ -278,7 +301,7 @@ struct heptaquark {
           if (std::abs(hqtrackd3.hqId()) != 3122)
             continue;
 
-          if (hqtrackd3.hqMass() < minLambdaMass || hqtrackd3.hqMass() > maxLambdaMass)
+          if (hqtrackd3.hqMass() < cfgMinLambdaMass || hqtrackd3.hqMass() > cfgMaxLambdaMass)
             continue;
 
           int isLambda = static_cast<int>(hqtrackd3.hqId() < 0);
@@ -321,6 +344,9 @@ struct heptaquark {
           HQ12 = HQ1 + HQ2;
           HQ13 = HQ1 + HQ3;
 
+          if (cfgSelectHQ && selectHQ(HQ1, HQ2, HQ3))
+            continue;
+
           histos.fill(HIST("h_InvMass_same"), exotic.M(), exotic.Pt(), collision.centrality());
           histos.fill(HIST("hDalitz"), HQ12.M2(), HQ13.M2(), exotic.M(), exotic.Pt(), isLambda, collision.centrality());