diff --git a/PWGHF/HFC/TableProducer/derivedDataCreatorCorrelationsReduced.cxx b/PWGHF/HFC/TableProducer/derivedDataCreatorCorrelationsReduced.cxx
index 3efa36a6d9f..73157c06d08 100644
--- a/PWGHF/HFC/TableProducer/derivedDataCreatorCorrelationsReduced.cxx
+++ b/PWGHF/HFC/TableProducer/derivedDataCreatorCorrelationsReduced.cxx
@@ -17,7 +17,6 @@
 
 #include "PWGHF/Core/CentralityEstimation.h"
 #include "PWGHF/Core/DecayChannels.h"
-#include "PWGHF/Core/HfHelper.h"
 #include "PWGHF/DataModel/AliasTables.h"
 #include "PWGHF/DataModel/CandidateReconstructionTables.h"
 #include "PWGHF/DataModel/CandidateSelectionTables.h"
@@ -33,7 +32,6 @@
 
 #include <CCDB/BasicCCDBManager.h>
 #include <CommonConstants/MathConstants.h>
-#include <CommonConstants/PhysicsConstants.h>
 #include <Framework/ASoA.h>
 #include <Framework/AnalysisDataModel.h>
 #include <Framework/AnalysisHelpers.h>
@@ -157,7 +155,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
   void init(InitContext&)
   {
 
-    hfEvSel.addHistograms(registry); // collision monitoring
+    hfEvSel.init(registry);
+
     ccdb->setURL(ccdbUrl);
     ccdb->setCaching(true);
     ccdb->setLocalObjectValidityChecking();
@@ -176,7 +175,7 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
 
-    hfcFlowHelper.init(ccdb, registry, classesMl, yCandMax.value);
+    hfcFlowHelper.init(classesMl, yCandMax.value);
 
     const AxisSpec axisTpcSignal = {binsTpcSignal, "TPC Signal"};
     const AxisSpec axisCent = {binsCent, "Centrality"};
@@ -209,32 +208,6 @@ struct HfDerivedDataCreatorCorrelationsReduced {
     }
   }; // end init
 
-  /// Check event selections for collision and fill the collision table
-  /// \param collision is the collision
-  template <typename Coll>
-  bool checkCollision(Coll const& collision, float& cent, float& mult)
-  {
-    o2::hf_evsel::HfCollisionRejectionMask collRejMask{};
-    if (centEstimator == CentralityEstimator::FT0A) {
-      collRejMask = hfEvSel.getHfCollisionRejectionMask<true, CentralityEstimator::FT0A, aod::BCsWithTimestamps>(collision, cent, ccdb, registry);
-      mult = collision.multFT0A();
-    } else if (centEstimator == CentralityEstimator::FT0C) {
-      collRejMask = hfEvSel.getHfCollisionRejectionMask<true, CentralityEstimator::FT0C, aod::BCsWithTimestamps>(collision, cent, ccdb, registry);
-      mult = collision.multFT0C();
-    } else if (centEstimator == CentralityEstimator::FT0M) {
-      collRejMask = hfEvSel.getHfCollisionRejectionMask<true, CentralityEstimator::FT0M, aod::BCsWithTimestamps>(collision, cent, ccdb, registry);
-      mult = collision.multFT0M();
-    } else if (centEstimator == CentralityEstimator::FV0A) {
-      collRejMask = hfEvSel.getHfCollisionRejectionMask<true, CentralityEstimator::FV0A, aod::BCsWithTimestamps>(collision, cent, ccdb, registry);
-      mult = collision.multFV0A();
-    } else {
-      LOG(fatal) << "Centrality estimator not recognized for collision selection";
-      std::abort();
-    }
-    hfEvSel.fillHistograms(collision, collRejMask, cent);
-    return collRejMask == 0;
-  }
-
   /// Checks if the trigger cand-associated track pair can be accepted for SE correlation
   /// \param assTrk is the associated track
   /// \param cand is the trigger candidate
@@ -404,7 +377,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -421,7 +395,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -441,7 +416,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -461,7 +437,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -482,7 +459,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -502,7 +480,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -523,7 +502,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -543,7 +523,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -564,7 +545,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -582,7 +564,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -602,7 +585,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -622,7 +606,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -643,7 +628,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -663,7 +649,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -684,7 +671,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -704,7 +692,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
       return;
     }
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -720,7 +709,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
                                  TracksData const& tracks)
   {
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -733,7 +723,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
                                  TracksData const& tracks)
   {
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -747,7 +738,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
                                TracksMc const& tracks)
   {
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -760,7 +752,8 @@ struct HfDerivedDataCreatorCorrelationsReduced {
                                TracksMc const& tracks)
   {
     float cent{-1.}, mult{-1.};
-    if (!checkCollision(coll, cent, mult)) {
+    auto centDetector = static_cast<CentralityEstimator>(centEstimator.value);
+    if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
       return;
     }
     rowCollisions(mult, coll.numContrib(), cent, coll.posZ());
@@ -769,7 +762,7 @@ struct HfDerivedDataCreatorCorrelationsReduced {
   PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcHadronHadronME, "Process Mixed Event for hadron candidates", false);
 };
 
-WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+WorkflowSpec defineDataProcessing(ConfigContext const& context)
 {
-  return WorkflowSpec{adaptAnalysisTask<HfDerivedDataCreatorCorrelationsReduced>(cfgc)};
+  return WorkflowSpec{adaptAnalysisTask<HfDerivedDataCreatorCorrelationsReduced>(context)};
 }
diff --git a/PWGHF/HFC/Tasks/taskCorrelationsEfficiencies.cxx b/PWGHF/HFC/Tasks/taskCorrelationsEfficiencies.cxx
index 8407dcbde82..365920f243b 100644
--- a/PWGHF/HFC/Tasks/taskCorrelationsEfficiencies.cxx
+++ b/PWGHF/HFC/Tasks/taskCorrelationsEfficiencies.cxx
@@ -14,13 +14,14 @@
 /// \author Marcello Di Costanzo <marcello.di.costanzo@cern.ch>, Politecnico and INFN Torino
 
 #include "PWGHF/Core/CentralityEstimation.h"
-#include "PWGHF/Core/HfHelper.h"
+#include "PWGHF/Core/DecayChannels.h"
 #include "PWGHF/DataModel/AliasTables.h"
 #include "PWGHF/DataModel/CandidateReconstructionTables.h"
 #include "PWGHF/DataModel/CandidateSelectionTables.h"
 #include "PWGHF/HFC/DataModel/DerivedDataCorrelationTables.h"
 #include "PWGHF/HFC/Utils/utilsFlow.h"
 #include "PWGHF/Utils/utilsEvSelHf.h"
+#include "PWGLF/DataModel/mcCentrality.h"
 
 #include "Common/Core/RecoDecay.h"
 #include "Common/DataModel/Centrality.h"
@@ -29,20 +30,22 @@
 #include "Common/DataModel/TrackSelectionTables.h"
 
 #include <CCDB/BasicCCDBManager.h>
-#include <CommonConstants/MathConstants.h>
-#include <CommonConstants/PhysicsConstants.h>
 #include <Framework/ASoA.h>
 #include <Framework/AnalysisDataModel.h>
 #include <Framework/AnalysisHelpers.h>
 #include <Framework/AnalysisTask.h>
 #include <Framework/Configurable.h>
+#include <Framework/DeviceSpec.h>
 #include <Framework/HistogramRegistry.h>
 #include <Framework/HistogramSpec.h>
 #include <Framework/InitContext.h>
 #include <Framework/Logger.h>
+#include <Framework/RunningWorkflowInfo.h>
 #include <Framework/runDataProcessing.h>
 
 #include <TF1.h>
+#include <THnSparse.h>
+#include <TPDGCode.h>
 #include <TString.h>
 
 #include <cmath>
@@ -127,13 +130,17 @@ struct HfTaskCorrelationsEfficiencies {
   Configurable<float> downSampleCharm{"downSampleCharm", 1., "Fraction of pairs to keep"};
   Configurable<float> ptMaxDownSampleCharm{"ptMaxDownSampleCharm", 10., "Maximum pt for the application of the downsampling factor"};
 
+  HfEventSelection hfEvSel; // event selection and monitoring
+  HfEventSelectionMc hfEvSelMc;
   HfcFlowHelper hfcFlowHelper; // flow utilities
   o2::framework::Service<o2::ccdb::BasicCCDBManager> ccdb{};
   SliceCache cache;
 
   std::unique_ptr<TF1> funcDcaXYPtCutPrimTrk{nullptr};
 
-  using RecoCollisions = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::PVMults>;
+  using BCsInfo = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels>;
+  using McCollisionsWMult = soa::Join<aod::McCollisions, aod::McCentFT0Ms>;
+  using RecoCollisions = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::PVMults, aod::FT0Mults, aod::FV0Mults, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFV0As>;
   using McCandsDs = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelDsToKKPi, aod::HfMlDsToKKPi, aod::HfCand3ProngMcRec>>;
   using McCandsDplus = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelDplusToPiKPi, aod::HfMlDplusToPiKPi, aod::HfCand3ProngMcRec>>;
   using McCandsD0 = soa::Filtered<soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfMlD0, aod::HfCand2ProngMcRec>>;
@@ -152,6 +159,7 @@ struct HfTaskCorrelationsEfficiencies {
   Preslice<McCandsLc> candsLcPerColl = aod::hf_cand::collisionId;
   Preslice<McTracks> trackIndicesPerColl = aod::track::collisionId;
 
+  PresliceUnsorted<aod::McParticles> mcPartsPerMcCollision = aod::mcparticle::mcCollisionId;
   PresliceUnsorted<aod::McCollisionLabels> colPerMcCollision = aod::mccollisionlabel::mcCollisionId;
 
   ConfigurableAxis binsPtCharm{"binsPtCharm", {VARIABLE_WIDTH, 0., 2., 4., 6., 8., 24.}, ""};
@@ -165,9 +173,18 @@ struct HfTaskCorrelationsEfficiencies {
 
   HistogramRegistry registry{"registry", {}};
 
-  void init(InitContext&)
+  void init(InitContext& initContext)
   {
 
+    hfEvSel.init(registry);
+    const auto& workflows = initContext.services().get<RunningWorkflowInfo const>();
+    for (const DeviceSpec& device : workflows.devices) {
+      if (device.name == "hf-task-correlations-efficiencies") {
+        // init HF event selection helper
+        hfEvSelMc.init(device, registry);
+        break;
+      }
+    }
     ccdb->setURL(ccdbUrl);
     ccdb->setCaching(true);
     ccdb->setLocalObjectValidityChecking();
@@ -186,7 +203,7 @@ struct HfTaskCorrelationsEfficiencies {
       return;
     }
 
-    hfcFlowHelper.init(ccdb, registry, classesMl, yCandMax.value);
+    hfcFlowHelper.init(classesMl, yCandMax.value);
 
     // Setup pt-dependent DCAxy cut function
     if (usePtDiffDcaXYCut) {
@@ -304,7 +321,8 @@ struct HfTaskCorrelationsEfficiencies {
         continue;
       }
       float cent{-1.}, mult{-1.};
-      if (!hfcFlowHelper.checkCollision<o2::hf_centrality::CentralityEstimator::None>(coll, cent, mult)) {
+      CentralityEstimator centDetector = CentralityEstimator::None;
+      if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
         continue;
       }
       fillCharmEfficiencies<TriggerChannels::DplusToPiKPi>(candsInColl, coll.numContrib());
@@ -357,20 +375,27 @@ struct HfTaskCorrelationsEfficiencies {
   void processEffAssTracks(RecoCollisions const& recoCollisions,
                            McTracks const& recoTracks,
                            aod::McParticles const& genParticles,
-                           aod::McCollisions const&)
+                           McCollisionsWMult const& genCollisions,
+                           BCsInfo const&)
   {
     for (const auto& coll : recoCollisions) {
       float cent{-1.}, mult{-1.};
-      if (!hfcFlowHelper.checkCollision<o2::hf_centrality::CentralityEstimator::None>(coll, cent, mult)) {
+      CentralityEstimator centDetector = CentralityEstimator::None;
+      if (!hfcFlowHelper.selectRecoCollision(hfEvSel, ccdb, registry, coll, centDetector, cent, mult)) {
         continue;
       }
       if (!coll.has_mcCollision()) {
         continue;
       }
-      auto mcCollision = coll.template mcCollision_as<aod::McCollisions>();
+      auto mcCollision = coll.template mcCollision_as<McCollisionsWMult>();
 
       auto tracksInColl = recoTracks.sliceBy(trackIndicesPerColl, coll.globalIndex());
       for (const auto& assTrk : tracksInColl) {
+
+        if (!assTrk.isGlobalTrackWoDCA() || assTrk.tpcNClsCrossedRows() < tpcNClsCrossedRowsMin) {
+          continue;
+        }
+
         if (!assTrk.has_mcParticle()) {
           continue;
         }
@@ -410,29 +435,39 @@ struct HfTaskCorrelationsEfficiencies {
     }
 
     // Fill histograms for generated associated tracks
-    for (const auto& genTrk : genParticles) {
-      int matchedSpecie = getAssocTrkSpecie(genTrk);
-      if (matchedSpecie < 0) {
-        continue;
-      }
-      if (!genTrk.has_mcCollision()) {
-        continue;
-      }
+    for (const auto& mcColl : genCollisions) {
 
-      // Pt and eta selection
-      if (genTrk.pt() < ptTrkMin || genTrk.pt() > ptTrkMax || std::fabs(genTrk.eta()) > etaTrkMax) {
+      float cent{-1.}, mult{-1.};
+      const auto recoColls = recoCollisions.sliceBy(colPerMcCollision, mcColl.globalIndex());
+      CentralityEstimator centDetector = CentralityEstimator::None;
+      if (!hfcFlowHelper.selectGenCollision<BCsInfo>(hfEvSelMc, mcColl, recoColls, centDetector, cent, mult)) {
         continue;
       }
 
-      int genTrkOrigin = genTrk.isPhysicalPrimary() ? 0 : 1; // 0 for primary, 1 for secondary
-      auto genCollision = genTrk.template mcCollision_as<aod::McCollisions>();
-      registry.fill(HIST("hAssTrackGen"), genTrk.pt(), genTrk.eta(), genTrk.phi(), genCollision.posZ(), matchedSpecie, genTrkOrigin);
+      const auto groupedMcParticles = genParticles.sliceBy(mcPartsPerMcCollision, mcColl.globalIndex());
+      for (const auto& genTrk : groupedMcParticles) {
+        int matchedSpecie = getAssocTrkSpecie(genTrk);
+        if (matchedSpecie < 0) {
+          continue;
+        }
+        if (!genTrk.has_mcCollision()) {
+          // Should never happen because we group on the generated collisions
+          continue;
+        }
+
+        // Pt and eta selection
+        if (genTrk.pt() < ptTrkMin || genTrk.pt() > ptTrkMax || std::fabs(genTrk.eta()) > etaTrkMax) {
+          continue;
+        }
+        int genTrkOrigin = genTrk.isPhysicalPrimary() ? 0 : 1; // 0 for primary, 1 for secondary
+        registry.fill(HIST("hAssTrackGen"), genTrk.pt(), genTrk.eta(), genTrk.phi(), mcColl.posZ(), matchedSpecie, genTrkOrigin);
+      }
     }
   }
   PROCESS_SWITCH(HfTaskCorrelationsEfficiencies, processEffAssTracks, "Process efficiencies for associated tracks", true);
 };
 
-WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+WorkflowSpec defineDataProcessing(ConfigContext const& context)
 {
-  return WorkflowSpec{adaptAnalysisTask<HfTaskCorrelationsEfficiencies>(cfgc)};
+  return WorkflowSpec{adaptAnalysisTask<HfTaskCorrelationsEfficiencies>(context)};
 }
diff --git a/PWGHF/HFC/Utils/utilsFlow.h b/PWGHF/HFC/Utils/utilsFlow.h
index 52b4885bfe4..8d345728769 100644
--- a/PWGHF/HFC/Utils/utilsFlow.h
+++ b/PWGHF/HFC/Utils/utilsFlow.h
@@ -17,31 +17,20 @@
 #define PWGHF_HFC_UTILS_UTILSFLOW_H_
 
 #include "PWGHF/Core/CentralityEstimation.h"
-#include "PWGHF/Core/DecayChannels.h"
 #include "PWGHF/Core/HfHelper.h"
 #include "PWGHF/Utils/utilsEvSelHf.h"
 
-#include "Common/Core/RecoDecay.h"
-
 #include <CCDB/BasicCCDBManager.h>
-#include <CommonConstants/MathConstants.h>
-#include <CommonConstants/PhysicsConstants.h>
 #include <Framework/ASoA.h>
 #include <Framework/AnalysisDataModel.h>
 #include <Framework/AnalysisHelpers.h>
-#include <Framework/AnalysisTask.h>
-#include <Framework/Configurable.h>
 #include <Framework/HistogramRegistry.h>
-#include <Framework/HistogramSpec.h>
-#include <Framework/InitContext.h>
 #include <Framework/Logger.h>
 #include <Framework/runDataProcessing.h>
 
-#include <TPDGCode.h>
-
 #include <cmath>
-#include <cstddef>
-#include <cstdint>
+#include <csignal>
+#include <cstdlib>
 #include <utility>
 #include <vector>
 
@@ -79,18 +68,12 @@ class HfcFlowHelper
   /// \param histRegistry is the histogram registry
   /// \param mlClasses are the machine learning classes
   /// \param maxRapidity is the maximum rapidity
-  void init(o2::framework::Service<o2::ccdb::BasicCCDBManager>& ccdbService,
-            o2::framework::HistogramRegistry& histRegistry,
-            std::vector<unsigned int> mlClasses,
+  void init(std::vector<unsigned int> mlClasses,
             float maxRapidity)
   {
-
-    LOG(info) << "Initializing HFC flow helper!";
-    this->hfEvSel.addHistograms(histRegistry);
-    this->ccdb = &ccdbService;
-    this->registry = &histRegistry;
     this->classMl = std::move(mlClasses);
     this->yCandMax = maxRapidity;
+    LOG(info) << "HFC flow helper initialized";
   }
 
   /// Get charm hadron candidate mass
@@ -168,6 +151,7 @@ class HfcFlowHelper
 
   /// Cut on rapidity of the candidate
   /// \param candidate is the charm hadron candidate
+  /// \return true if the candidate passes the rapidity cut, false otherwise
   template <TriggerChannels CandType, typename TCand>
   bool cutCandRapidity(const TCand& candidate)
   {
@@ -186,40 +170,103 @@ class HfcFlowHelper
     return std::fabs(y) < yCandMax;
   }
 
-  /// Check event selections for collision and fill the collision table
+  /// Runtime dispatch to the compile-time checkCollision
   /// \param collision is the collision
-  template <o2::hf_centrality::CentralityEstimator CentEst, typename Coll>
-  bool checkCollision(Coll const& collision, float& cent, float& mult)
+  /// \param centEst is the centrality estimator
+  /// \param cent is the centrality value
+  /// \param mult is the multiplicity value
+  template <typename BCsInfo, typename GenColl, typename RecoCollSlice>
+  bool selectGenCollision(o2::hf_evsel::HfEventSelectionMc& hfEvSelMc,
+                          GenColl const& genColl,
+                          RecoCollSlice const& recoCollsSlice,
+                          o2::hf_centrality::CentralityEstimator centEst,
+                          float& cent,
+                          float& mult)
   {
+    int nSplitColl = recoCollsSlice.size();
     o2::hf_evsel::HfCollisionRejectionMask collRejMask{};
-    if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::FT0A) {
-      collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::FT0A, aod::BCsWithTimestamps>(collision, cent, *ccdb, *registry);
-      mult = collision.multFT0A();
-    } else if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::FT0C) {
-      collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::FT0C, aod::BCsWithTimestamps>(collision, cent, *ccdb, *registry);
-      mult = collision.multFT0C();
-    } else if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::FT0M) {
-      collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::FT0M, aod::BCsWithTimestamps>(collision, cent, *ccdb, *registry);
-      mult = collision.multFT0M();
-    } else if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::FV0A) {
-      collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::FV0A, aod::BCsWithTimestamps>(collision, cent, *ccdb, *registry);
-      mult = collision.multFV0A();
-    } else if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::None) {
-      collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::None, aod::BCsWithTimestamps>(collision, cent, *ccdb, *registry);
-      mult = -1.f;
-    } else {
-      LOG(fatal) << "Centrality estimator not recognized for collision selection";
-      std::abort();
+    switch (centEst) {
+      case o2::hf_centrality::CentralityEstimator::FT0A:
+        collRejMask = hfEvSelMc.getHfMcCollisionRejectionMask<BCsInfo, o2::hf_centrality::CentralityEstimator::FT0A>(genColl, recoCollsSlice, cent);
+        mult = o2::hf_centrality::getCentralityGenColl(recoCollsSlice, o2::hf_centrality::CentralityEstimator::FT0A);
+        hfEvSelMc.fillHistograms<o2::hf_centrality::CentralityEstimator::FT0A>(genColl, collRejMask, nSplitColl);
+        break;
+      case o2::hf_centrality::CentralityEstimator::FT0C:
+        collRejMask = hfEvSelMc.getHfMcCollisionRejectionMask<BCsInfo, o2::hf_centrality::CentralityEstimator::FT0C>(genColl, recoCollsSlice, cent);
+        mult = o2::hf_centrality::getCentralityGenColl(recoCollsSlice, o2::hf_centrality::CentralityEstimator::FT0C);
+        hfEvSelMc.fillHistograms<o2::hf_centrality::CentralityEstimator::FT0C>(genColl, collRejMask, nSplitColl);
+        break;
+      case o2::hf_centrality::CentralityEstimator::FT0M:
+        collRejMask = hfEvSelMc.getHfMcCollisionRejectionMask<BCsInfo, o2::hf_centrality::CentralityEstimator::FT0M>(genColl, recoCollsSlice, cent);
+        mult = o2::hf_centrality::getCentralityGenColl(recoCollsSlice, o2::hf_centrality::CentralityEstimator::FT0M);
+        hfEvSelMc.fillHistograms<o2::hf_centrality::CentralityEstimator::FT0M>(genColl, collRejMask, nSplitColl);
+        break;
+      case o2::hf_centrality::CentralityEstimator::FV0A:
+        collRejMask = hfEvSelMc.getHfMcCollisionRejectionMask<BCsInfo, o2::hf_centrality::CentralityEstimator::FV0A>(genColl, recoCollsSlice, cent);
+        mult = o2::hf_centrality::getCentralityGenColl(recoCollsSlice, o2::hf_centrality::CentralityEstimator::FV0A);
+        hfEvSelMc.fillHistograms<o2::hf_centrality::CentralityEstimator::FV0A>(genColl, collRejMask, nSplitColl);
+        break;
+      case o2::hf_centrality::CentralityEstimator::None:
+        collRejMask = hfEvSelMc.getHfMcCollisionRejectionMask<BCsInfo, o2::hf_centrality::CentralityEstimator::None>(genColl, recoCollsSlice, cent);
+        mult = -1.f;
+        hfEvSelMc.fillHistograms<o2::hf_centrality::CentralityEstimator::None>(genColl, collRejMask, nSplitColl);
+        break;
+      default:
+        LOG(fatal) << "Centrality estimator not recognized for collision selection";
+        std::abort();
     }
+
+    return (collRejMask == 0);
+  }
+
+  /// Apply event selection on reconstructed collisions
+  /// \param collision is the collision
+  /// \param centEst is the centrality estimator
+  /// \param cent is the centrality value
+  /// \param mult is the multiplicity value
+  template <typename Coll>
+  bool selectRecoCollision(o2::hf_evsel::HfEventSelection& hfEvSel,
+                           o2::framework::Service<o2::ccdb::BasicCCDBManager>& ccdb,
+                           o2::framework::HistogramRegistry& registry,
+                           Coll const& collision,
+                           o2::hf_centrality::CentralityEstimator centEst,
+                           float& cent,
+                           float& mult)
+  {
+    o2::hf_evsel::HfCollisionRejectionMask collRejMask{};
+    switch (centEst) {
+      case o2::hf_centrality::CentralityEstimator::FT0A:
+        collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::FT0A, aod::BCsWithTimestamps>(collision, cent, ccdb, registry);
+        mult = collision.multFT0A();
+        break;
+      case o2::hf_centrality::CentralityEstimator::FT0C:
+        collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::FT0C, aod::BCsWithTimestamps>(collision, cent, ccdb, registry);
+        mult = collision.multFT0C();
+        break;
+      case o2::hf_centrality::CentralityEstimator::FT0M:
+        collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::FT0M, aod::BCsWithTimestamps>(collision, cent, ccdb, registry);
+        mult = collision.multFT0M();
+        break;
+      case o2::hf_centrality::CentralityEstimator::FV0A:
+        collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::FV0A, aod::BCsWithTimestamps>(collision, cent, ccdb, registry);
+        mult = collision.multFV0A();
+        break;
+      case o2::hf_centrality::CentralityEstimator::None:
+        collRejMask = hfEvSel.getHfCollisionRejectionMask<true, o2::hf_centrality::CentralityEstimator::None, aod::BCsWithTimestamps>(collision, cent, ccdb, registry);
+        mult = -1.f;
+        break;
+      default:
+        LOG(fatal) << "Centrality estimator not recognized for collision selection";
+        std::abort();
+    }
+
     hfEvSel.fillHistograms(collision, collRejMask, cent);
-    return collRejMask == 0;
+
+    return (collRejMask == 0);
   }
 
  private:
-  o2::hf_evsel::HfEventSelection hfEvSel;
-  o2::framework::Service<o2::ccdb::BasicCCDBManager>* ccdb{nullptr};
-  o2::framework::HistogramRegistry* registry{nullptr};
-  std::vector<unsigned int> classMl{};
+  std::vector<unsigned int> classMl;
   float yCandMax{-999.f};
 };