diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx
index 14d9c2c190f..6e228d572dc 100644
--- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx
+++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx
@@ -105,6 +105,12 @@ enum DecayChannel { DplusToPiKPi = 0,
                     Xic0ToXiPi
 };
 
+enum RunMode {
+  kSP = 0,
+  kEsE,
+  kSPEsE,
+};
+
 struct HfTaskFlowCharmHadrons {
   Produces<o2::aod::HfCandMPtInfos> rowCandMassPtMl;
   Produces<o2::aod::HfCandFlowInfos> rowCandMassPtMlSpCent;
@@ -119,6 +125,7 @@ struct HfTaskFlowCharmHadrons {
   Configurable<float> centralityMax{"centralityMax", 100., "Maximum centrality accepted in SP/EP computation (not applied in resolution process)"};
   Configurable<bool> storeEP{"storeEP", false, "Flag to store EP-related axis"};
   Configurable<bool> storeMl{"storeMl", false, "Flag to store ML scores"};
+  Configurable<bool> storeSPQVec{"storeSPQVec", true, "Flag to store the Q-vectors for SP"};
   Configurable<bool> storeRedQVec{"storeRedQVec", false, "Flag to store reduced Q-vectors for ESE"};
   Configurable<bool> fillMassPtMlTree{"fillMassPtMlTree", false, "Flag to fill mass, pt and ML scores tree"};
   Configurable<bool> fillMassPtMlSpCentTree{"fillMassPtMlSpCentTree", false, "Flag to fill mass, pt, ML scores, SP and centrality tree"};
@@ -153,7 +160,9 @@ struct HfTaskFlowCharmHadrons {
   using CandXic0DataWMl = soa::Filtered<soa::Join<aod::HfCandToXiPiKf, aod::HfSelToXiPiKf, aod::HfMlToXiPi>>;
   using CandD0DataWMl = soa::Filtered<soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfMlD0>>;
   using CandD0Data = soa::Filtered<soa::Join<aod::HfCand2Prong, aod::HfSelD0>>;
-  using CollsWithQvecs = soa::Join<aod::Collisions, aod::EvSels, aod::QvectorFT0Cs, aod::QvectorFT0As, aod::QvectorFT0Ms, aod::QvectorFV0As, aod::QvectorTPCposs, aod::QvectorTPCnegs, aod::QvectorTPCalls, aod::EseQvecPercs, aod::CentFV0As, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0CVariant2s>;
+  using CollsWithSPQvecs = soa::Join<aod::Collisions, aod::EvSels, aod::QvectorFT0Cs, aod::QvectorFT0As, aod::QvectorFT0Ms, aod::QvectorFV0As, aod::QvectorBPoss, aod::QvectorBNegs, aod::QvectorBTots, aod::CentFV0As, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0CVariant2s>;
+  using CollsWithEsEQvecs = soa::Join<aod::Collisions, aod::EvSels, aod::EseQvecPercs, aod::CentFV0As, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0CVariant2s>;
+  using CollsWithSPEsEQvecs = soa::Join<aod::Collisions, aod::EvSels, aod::EseQvecPercs, aod::QvectorFT0Cs, aod::QvectorFT0As, aod::QvectorFT0Ms, aod::QvectorFV0As, aod::QvectorBPoss, aod::QvectorBNegs, aod::QvectorBTots, aod::CentFV0As, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0CVariant2s>;
   using TracksWithExtra = soa::Join<aod::Tracks, aod::TracksExtra>;
 
   Filter filterSelectDsCandidates = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag || aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag;
@@ -212,6 +221,12 @@ struct HfTaskFlowCharmHadrons {
     if (storeResoOccu && occEstimator == 0) {
       LOGP(fatal, "Occupancy estimation must be enabled to store resolution THnSparse! Please check your configuration!");
     }
+    if (storeRedQVec && !doprocessD0EsEMl && !doprocessD0SPEsEMl) {
+      LOGP(fatal, "EsE q-vectors require dedicated process functions. Please check your configuration!");
+    }
+    if (storeSPQVec && doprocessD0EsEMl) {
+      LOGP(fatal, "Scalar product tables not available for this process function. Please check your configuration!");
+    }
     const AxisSpec thnAxisInvMass{thnConfigAxisInvMass, "Inv. mass (GeV/#it{c}^{2})"};
     const AxisSpec thnAxisPt{thnConfigAxisPt, "#it{p}_{T} (GeV/#it{c})"};
     const AxisSpec thnAxisCent{thnConfigAxisCent, "Centrality"};
@@ -237,7 +252,10 @@ struct HfTaskFlowCharmHadrons {
     const AxisSpec thnAxisResoFT0cTPCtot{thnConfigAxisResoFT0cTPCtot, "Q_{FT0c} #bullet Q_{TPCtot}"};
     const AxisSpec thnAxisResoFV0aTPCtot{thnConfigAxisResoFV0aTPCtot, "Q_{FV0a} #bullet Q_{TPCtot}"};
 
-    std::vector<AxisSpec> axes = {thnAxisInvMass, thnAxisPt, thnAxisCent, thnAxisScalarProd};
+    std::vector<AxisSpec> axes = {thnAxisInvMass, thnAxisPt, thnAxisCent};
+    if (storeSPQVec) {
+      axes.insert(axes.end(), {thnAxisScalarProd});
+    }
     if (storeEP) {
       axes.insert(axes.end(), {thnAxisCosNPhi, thnAxisSinNPhi, thnAxisCosDeltaPhi});
     }
@@ -279,7 +297,7 @@ struct HfTaskFlowCharmHadrons {
       registry.add("ep/hSparseEp", "THn for Event Plane distirbution", {HistType::kTHnSparseF, {thnAxisCent, thnAxisPsi, thnAxisCosNPhi, thnAxisSinNPhi}});
     }
 
-    if (doprocessResolution) { // enable resolution histograms only for resolution process
+    if (doprocessResolutionSP || doprocessResolutionSPEsE) { // enable resolution histograms only for resolution process
       registry.add("spReso/hSpResoFT0cFT0a", "hSpResoFT0cFT0a; centrality; Q_{FT0c} #bullet Q_{FT0a}", {HistType::kTH2F, {thnAxisCent, thnAxisScalarProd}});
       registry.add("spReso/hSpResoFT0cFV0a", "hSpResoFT0cFV0a; centrality; Q_{FT0c} #bullet Q_{FV0a}", {HistType::kTH2F, {thnAxisCent, thnAxisScalarProd}});
       registry.add("spReso/hSpResoFT0cTPCpos", "hSpResoFT0cTPCpos; centrality; Q_{FT0c} #bullet Q_{TPCpos}", {HistType::kTH2F, {thnAxisCent, thnAxisScalarProd}});
@@ -298,9 +316,12 @@ struct HfTaskFlowCharmHadrons {
       registry.add("spReso/hSpResoFV0aTPCtot", "hSpResoFV0aTPCtot; centrality; Q_{FV0a} #bullet Q_{TPCtot}", {HistType::kTH2F, {thnAxisCent, thnAxisScalarProd}});
       registry.add("spReso/hSpResoTPCposTPCneg", "hSpResoTPCposTPCneg; centrality; Q_{TPCpos} #bullet Q_{TPCneg}", {HistType::kTH2F, {thnAxisCent, thnAxisScalarProd}});
 
-      if (storeRedQVec) {
+      if (doprocessResolutionSPEsE) {
         registry.add("redQVecs/hSparseRedQVecs", "hSpResoRedQVec; centrality; Q_{red} #bullet Q_{TPCtot}", {HistType::kTHnSparseF, {thnAxisCent, thnAxisRedQVec, thnAxisRedQVec, thnAxisRedQVec, thnAxisRedQVec}});
         registry.add("redQVecs/hRedQVecFT0C", "hRedQVecFT0C; centrality; q_{FT0c}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}});
+        registry.add("redQVecs/hRedQVecFT0A", "hRedQVecFT0A; centrality; q_{FT0a}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}});
+        registry.add("redQVecs/hRedQVecFT0M", "hRedQVecFT0M; centrality; q_{FT0m}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}});
+        registry.add("redQVecs/hRedQVecFV0A", "hRedQVecFV0A; centrality; q_{FV0a}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}});
         registry.add("redQVecs/hRedQVecTpcPos", "hRedQVecTpcPos; centrality; q_{TPCpos}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}});
         registry.add("redQVecs/hRedQVecTpcNeg", "hRedQVecTpcNeg; centrality; q_{TPCneg}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}});
         registry.add("redQVecs/hRedQVecTpcAll", "hRedQVecTpcAll; centrality; q_{TPCall}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}});
@@ -472,7 +493,9 @@ struct HfTaskFlowCharmHadrons {
     values.push_back(mass);
     values.push_back(pt);
     values.push_back(cent);
-    values.push_back(sp);
+    if (storeSPQVec) {
+      values.push_back(sp);
+    }
 
     if (storeEP) {
       values.push_back(cosNPhi);
@@ -518,8 +541,8 @@ struct HfTaskFlowCharmHadrons {
   /// \param bc is the bunch crossing with timestamp information
   /// \param centrality is the collision centrality
   /// \return true if the collision is selected, false otherwise
-  template <o2::hf_centrality::CentralityEstimator CentEstimator>
-  bool isCollSelected(CollsWithQvecs::iterator const& collision,
+  template <o2::hf_centrality::CentralityEstimator CentEstimator, typename TCollision>
+  bool isCollSelected(TCollision const& collision,
                       aod::BCsWithTimestamps const&,
                       float& centrality)
   {
@@ -539,8 +562,8 @@ struct HfTaskFlowCharmHadrons {
   /// Compute the scalar product
   /// \param collision is the collision with the Q vector information and event plane
   /// \param candidates are the selected candidates
-  template <DecayChannel Channel, typename T1, typename Trk>
-  void runFlowAnalysis(CollsWithQvecs::iterator const& collision,
+  template <RunMode StoreInfo, DecayChannel Channel, typename TCollision, typename T1, typename Trk>
+  void runFlowAnalysis(TCollision const& collision,
                        T1 const& candidates,
                        Trk const& /*tracks*/)
   {
@@ -563,7 +586,10 @@ struct HfTaskFlowCharmHadrons {
     float const sigmaMDplus = 0.015; // used 15 MeV as the D+ average peak width in run3
     float const nSigmaMass = 2.5;
 
-    std::array<float, 3> qVecs = getQvec(collision, qVecDetector.value);
+    std::array<float, 3> qVecs{-999.f, -999.f, -999.f};
+    if constexpr (StoreInfo == RunMode::kSP || StoreInfo == RunMode::kSPEsE) {
+      qVecs = getQvec(collision, qVecDetector.value);
+    }
     float xQVec = qVecs[0];
     float yQVec = qVecs[1];
     float const amplQVec = qVecs[2];
@@ -571,7 +597,7 @@ struct HfTaskFlowCharmHadrons {
 
     float redQVec{-999.f};
     std::array<float, 3> qVecRedComps{-999.f, -999.f, -999.f};
-    if (storeRedQVec) {
+    if constexpr (StoreInfo == RunMode::kEsE || StoreInfo == RunMode::kSPEsE) {
       qVecRedComps = getEseQvec(collision, qVecRedDetector.value);
     }
     float xRedQVec = qVecRedComps[0];
@@ -798,142 +824,169 @@ struct HfTaskFlowCharmHadrons {
   }
 
   // Ds with ML
-  void processDsMl(CollsWithQvecs::iterator const& collision,
+  void processDsMl(CollsWithSPQvecs::iterator const& collision,
                    CandDsDataWMl const&,
                    TracksWithExtra const& tracks)
   {
     auto candsDsToKKPiWMl = selectedDsToKKPiWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
     auto candsDsToPiKKWMl = selectedDsToPiKKWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::DsToKKPi>(collision, candsDsToKKPiWMl, tracks);
-    runFlowAnalysis<DecayChannel::DsToPiKK>(collision, candsDsToPiKKWMl, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::DsToKKPi>(collision, candsDsToKKPiWMl, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::DsToPiKK>(collision, candsDsToPiKKWMl, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processDsMl, "Process Ds candidates with ML", false);
 
   // Ds with rectangular cuts
-  void processDs(CollsWithQvecs::iterator const& collision,
+  void processDs(CollsWithSPQvecs::iterator const& collision,
                  CandDsData const& /*candidatesDs*/,
                  TracksWithExtra const& tracks)
   {
     auto candsDsToKKPi = selectedDsToKKPi->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
     auto candsDsToPiKK = selectedDsToPiKK->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::DsToKKPi>(collision, candsDsToKKPi, tracks);
-    runFlowAnalysis<DecayChannel::DsToPiKK>(collision, candsDsToPiKK, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::DsToKKPi>(collision, candsDsToKKPi, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::DsToPiKK>(collision, candsDsToPiKK, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processDs, "Process Ds candidates", false);
 
   // Dplus with ML
-  void processDplusMl(CollsWithQvecs::iterator const& collision,
+  void processDplusMl(CollsWithSPQvecs::iterator const& collision,
                       CandDplusDataWMl const& candidatesDplus,
                       TracksWithExtra const& tracks)
   {
-    runFlowAnalysis<DecayChannel::DplusToPiKPi>(collision, candidatesDplus, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::DplusToPiKPi>(collision, candidatesDplus, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processDplusMl, "Process Dplus candidates with ML", false);
 
   // Dplus with rectangular cuts
-  void processDplus(CollsWithQvecs::iterator const& collision,
+  void processDplus(CollsWithSPQvecs::iterator const& collision,
                     CandDplusData const& candidatesDplus,
                     TracksWithExtra const& tracks)
   {
-    runFlowAnalysis<DecayChannel::DplusToPiKPi>(collision, candidatesDplus, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::DplusToPiKPi>(collision, candidatesDplus, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processDplus, "Process Dplus candidates", true);
 
   // D0 with ML
-  void processD0Ml(CollsWithQvecs::iterator const& collision,
+  void processD0Ml(CollsWithSPQvecs::iterator const& collision,
                    CandD0DataWMl const& /*candidatesD0*/,
                    TracksWithExtra const& tracks)
   {
     auto candsD0ToPiKWMl = selectedD0ToPiKWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
     auto candsD0ToKPiWMl = selectedD0ToKPiWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::D0ToPiK>(collision, candsD0ToPiKWMl, tracks);
-    runFlowAnalysis<DecayChannel::D0ToKPi>(collision, candsD0ToKPiWMl, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::D0ToPiK>(collision, candsD0ToPiKWMl, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::D0ToKPi>(collision, candsD0ToKPiWMl, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processD0Ml, "Process D0 candidates with ML", false);
 
+  // D0 with ML
+  void processD0EsEMl(CollsWithEsEQvecs::iterator const& collision,
+                      CandD0DataWMl const& /*candidatesD0*/,
+                      TracksWithExtra const& tracks)
+  {
+    auto candsD0ToPiKWMl = selectedD0ToPiKWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
+    auto candsD0ToKPiWMl = selectedD0ToKPiWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
+    runFlowAnalysis<RunMode::kEsE, DecayChannel::D0ToPiK>(collision, candsD0ToPiKWMl, tracks);
+    runFlowAnalysis<RunMode::kEsE, DecayChannel::D0ToKPi>(collision, candsD0ToKPiWMl, tracks);
+  }
+  PROCESS_SWITCH(HfTaskFlowCharmHadrons, processD0EsEMl, "Process D0 candidates with ML for EsE information", false);
+
+  // D0 with ML
+  void processD0SPEsEMl(CollsWithSPEsEQvecs::iterator const& collision,
+                        CandD0DataWMl const& /*candidatesD0*/,
+                        TracksWithExtra const& tracks)
+  {
+    auto candsD0ToPiKWMl = selectedD0ToPiKWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
+    auto candsD0ToKPiWMl = selectedD0ToKPiWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
+    runFlowAnalysis<RunMode::kSPEsE, DecayChannel::D0ToPiK>(collision, candsD0ToPiKWMl, tracks);
+    runFlowAnalysis<RunMode::kSPEsE, DecayChannel::D0ToKPi>(collision, candsD0ToKPiWMl, tracks);
+  }
+  PROCESS_SWITCH(HfTaskFlowCharmHadrons, processD0SPEsEMl, "Process D0 candidates with ML for SP and EsE information", false);
+
   // D0 with rectangular cuts
-  void processD0(CollsWithQvecs::iterator const& collision,
+  void processD0(CollsWithSPQvecs::iterator const& collision,
                  CandD0Data const& /*candidatesD0*/,
                  TracksWithExtra const& tracks)
   {
     auto candsD0ToPiK = selectedD0ToPiK->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
     auto candsD0ToKPi = selectedD0ToKPi->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::D0ToPiK>(collision, candsD0ToPiK, tracks);
-    runFlowAnalysis<DecayChannel::D0ToKPi>(collision, candsD0ToKPi, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::D0ToPiK>(collision, candsD0ToPiK, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::D0ToKPi>(collision, candsD0ToKPi, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processD0, "Process D0 candidates", false);
 
   // Lc with ML
-  void processLcMl(CollsWithQvecs::iterator const& collision,
+  void processLcMl(CollsWithSPQvecs::iterator const& collision,
                    CandLcDataWMl const& /*candidatesLc*/,
                    TracksWithExtra const& tracks)
   {
     auto candsLcToPKPiWMl = selectedLcToPKPiWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
     auto candsLcToPiKPWMl = selectedLcToPiKPWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::LcToPKPi>(collision, candsLcToPKPiWMl, tracks);
-    runFlowAnalysis<DecayChannel::LcToPiKP>(collision, candsLcToPiKPWMl, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::LcToPKPi>(collision, candsLcToPKPiWMl, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::LcToPiKP>(collision, candsLcToPiKPWMl, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processLcMl, "Process Lc candidates with ML", false);
 
   // Lc with rectangular cuts
-  void processLc(CollsWithQvecs::iterator const& collision,
+  void processLc(CollsWithSPQvecs::iterator const& collision,
                  CandLcData const& /*candidatesLc*/,
                  TracksWithExtra const& tracks)
   {
     auto candsLcToPKPi = selectedLcToPKPi->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
     auto candsLcToPiKP = selectedLcToPiKP->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::LcToPKPi>(collision, candsLcToPKPi, tracks);
-    runFlowAnalysis<DecayChannel::LcToPiKP>(collision, candsLcToPiKP, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::LcToPKPi>(collision, candsLcToPKPi, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::LcToPiKP>(collision, candsLcToPiKP, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processLc, "Process Lc candidates", false);
 
   // Xic with ML
-  void processXicMl(CollsWithQvecs::iterator const& collision,
+  void processXicMl(CollsWithSPQvecs::iterator const& collision,
                     CandXicDataWMl const& /*candidatesXic*/,
                     TracksWithExtra const& tracks)
   {
     auto candsXicToPKPiWMl = selectedXicToPKPiWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
     auto candsXicToPiKPWMl = selectedXicToPiKPWMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::XicToPKPi>(collision, candsXicToPKPiWMl, tracks);
-    runFlowAnalysis<DecayChannel::XicToPiKP>(collision, candsXicToPiKPWMl, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::XicToPKPi>(collision, candsXicToPKPiWMl, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::XicToPiKP>(collision, candsXicToPiKPWMl, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processXicMl, "Process Xic candidates with ML", false);
 
   // Xic with rectangular cuts
-  void processXic(CollsWithQvecs::iterator const& collision,
+  void processXic(CollsWithSPQvecs::iterator const& collision,
                   CandXicData const& /*candidatesXic*/,
                   TracksWithExtra const& tracks)
   {
     auto candsXicToPKPi = selectedXicToPKPi->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
     auto candsXicToPiKP = selectedXicToPiKP->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::XicToPKPi>(collision, candsXicToPKPi, tracks);
-    runFlowAnalysis<DecayChannel::XicToPiKP>(collision, candsXicToPiKP, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::XicToPKPi>(collision, candsXicToPKPi, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::XicToPiKP>(collision, candsXicToPiKP, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processXic, "Process Xic candidates", false);
 
   // Xic0 with ML
-  void processXic0Ml(CollsWithQvecs::iterator const& collision,
+  void processXic0Ml(CollsWithSPQvecs::iterator const& collision,
                      CandXic0DataWMl const& /*candidatesXic0*/,
                      TracksWithExtra const& tracks)
   {
     auto candsXic0WMl = selectedXic0WMl->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::Xic0ToXiPi>(collision, candsXic0WMl, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::Xic0ToXiPi>(collision, candsXic0WMl, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processXic0Ml, "Process Xic0 candidates with ML", false);
 
   // Xic0
-  void processXic0(CollsWithQvecs::iterator const& collision,
+  void processXic0(CollsWithSPQvecs::iterator const& collision,
                    CandXic0Data const& /*candidatesXic0*/,
                    TracksWithExtra const& tracks)
   {
     auto candsXic0 = selectedXic0->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
-    runFlowAnalysis<DecayChannel::Xic0ToXiPi>(collision, candsXic0, tracks);
+    runFlowAnalysis<RunMode::kSP, DecayChannel::Xic0ToXiPi>(collision, candsXic0, tracks);
   }
   PROCESS_SWITCH(HfTaskFlowCharmHadrons, processXic0, "Process Xic0 candidates", false);
 
-  // Resolution
-  void processResolution(CollsWithQvecs::iterator const& collision,
-                         aod::BCsWithTimestamps const& bcs)
+  /// Compute resolution
+  /// \param collision is the collision with the Q vector information
+  /// \param bcs is used for the event selection based on centrality estimators with time information
+  template <bool HasRedQVecs, typename TCollision>
+  void fillResolution(TCollision const& collision,
+                      aod::BCsWithTimestamps const& bcs)
   {
     float centrality{-1.f};
     float const xQVecFT0a = collision.qvecFT0ARe();
@@ -944,12 +997,12 @@ struct HfTaskFlowCharmHadrons {
     float const yQVecFT0m = collision.qvecFT0MIm();
     float const xQVecFV0a = collision.qvecFV0ARe();
     float const yQVecFV0a = collision.qvecFV0AIm();
-    float const xQVecTPCPos = collision.qvecTPCposRe();
-    float const yQVecTPCPos = collision.qvecTPCposIm();
-    float const xQVecTPCNeg = collision.qvecTPCnegRe();
-    float const yQVecTPCNeg = collision.qvecTPCnegIm();
-    float const xQVecTPCAll = collision.qvecTPCallRe();
-    float const yQVecTPCAll = collision.qvecTPCallIm();
+    float const xQVecTPCPos = collision.qvecBPosRe();
+    float const yQVecTPCPos = collision.qvecBPosIm();
+    float const xQVecTPCNeg = collision.qvecBNegRe();
+    float const yQVecTPCNeg = collision.qvecBNegIm();
+    float const xQVecTPCAll = collision.qvecBTotRe();
+    float const yQVecTPCAll = collision.qvecBTotIm();
 
     centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator);
     if (storeResoOccu) {
@@ -970,7 +1023,7 @@ struct HfTaskFlowCharmHadrons {
       registry.fill(HIST("hSparseCentEstimators"), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(0)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(1)),
                     o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(2)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(3)));
     }
-    if (storeRedQVec) {
+    if constexpr (HasRedQVecs) {
       registry.fill(HIST("redQVecs/hSparseRedQVecs"), centrality,
                     std::hypot(collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm()),
                     std::hypot(collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm()),
@@ -1001,11 +1054,15 @@ struct HfTaskFlowCharmHadrons {
     registry.fill(HIST("spReso/hSpResoFV0aTPCtot"), centrality, xQVecFV0a * xQVecTPCAll + yQVecFV0a * yQVecTPCAll);
     registry.fill(HIST("spReso/hSpResoTPCposTPCneg"), centrality, xQVecTPCPos * xQVecTPCNeg + yQVecTPCPos * yQVecTPCNeg);
 
-    registry.fill(HIST("redQVecs/hRedQVecFT0C"), centrality, std::hypot(collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm()));
-    registry.fill(HIST("redQVecs/hRedQVecTpcPos"), centrality, std::hypot(collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm()));
-    registry.fill(HIST("redQVecs/hRedQVecTpcNeg"), centrality, std::hypot(collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm()));
-    registry.fill(HIST("redQVecs/hRedQVecTpcAll"), centrality, std::hypot(collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm()));
-
+    if constexpr (HasRedQVecs) {
+      registry.fill(HIST("redQVecs/hRedQVecFT0C"), centrality, std::hypot(collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm()));
+      registry.fill(HIST("redQVecs/hRedQVecFT0M"), centrality, std::hypot(collision.eseQvecFT0MRe(), collision.eseQvecFT0MIm()));
+      registry.fill(HIST("redQVecs/hRedQVecFT0A"), centrality, std::hypot(collision.eseQvecFT0ARe(), collision.eseQvecFT0AIm()));
+      registry.fill(HIST("redQVecs/hRedQVecFV0A"), centrality, std::hypot(collision.eseQvecFV0ARe(), collision.eseQvecFV0AIm()));
+      registry.fill(HIST("redQVecs/hRedQVecTpcPos"), centrality, std::hypot(collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm()));
+      registry.fill(HIST("redQVecs/hRedQVecTpcNeg"), centrality, std::hypot(collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm()));
+      registry.fill(HIST("redQVecs/hRedQVecTpcAll"), centrality, std::hypot(collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm()));
+    }
     if (saveEpResoHisto) {
       float const epFT0a = epHelper.GetEventPlane(xQVecFT0a, yQVecFT0a, harmonic);
       float const epFT0c = epHelper.GetEventPlane(xQVecFT0c, yQVecFT0c, harmonic);
@@ -1041,7 +1098,20 @@ struct HfTaskFlowCharmHadrons {
                     std::sin(harmonic * epHelper.GetEventPlane(xQVecFT0c, yQVecFT0c, harmonic)));
     }
   }
-  PROCESS_SWITCH(HfTaskFlowCharmHadrons, processResolution, "Process resolution", false);
+
+  void processResolutionSP(CollsWithSPQvecs::iterator const& collision,
+                           aod::BCsWithTimestamps const& bcs)
+  {
+    fillResolution<false>(collision, bcs);
+  }
+  PROCESS_SWITCH(HfTaskFlowCharmHadrons, processResolutionSP, "Process SP resolution", false);
+
+  void processResolutionSPEsE(CollsWithSPEsEQvecs::iterator const& collision,
+                              aod::BCsWithTimestamps const& bcs)
+  {
+    fillResolution<true>(collision, bcs);
+  }
+  PROCESS_SWITCH(HfTaskFlowCharmHadrons, processResolutionSPEsE, "Process SP and EsE resolution", false);
 
 }; // End struct HfTaskFlowCharmHadrons
 
diff --git a/PWGHF/D2H/Utils/utilsFlow.h b/PWGHF/D2H/Utils/utilsFlow.h
index 911fc60dd2c..02171203eba 100644
--- a/PWGHF/D2H/Utils/utilsFlow.h
+++ b/PWGHF/D2H/Utils/utilsFlow.h
@@ -202,6 +202,49 @@ std::array<float, 3> getQvec(TCollision const& collision, const int qvecEst)
   return std::array<float, 3>{-999.f, -999.f, -999.f};
 }
 
+/// Get the EsE Q vector
+template <typename T>
+concept HasEsEQvecFT0A = requires(T collision) {
+  collision.eseQvecFT0ARe();
+  collision.eseQvecFT0AIm();
+};
+
+template <typename T>
+concept HasEsEQvecFT0C = requires(T collision) {
+  collision.eseQvecFT0CRe();
+  collision.eseQvecFT0CIm();
+};
+
+template <typename T>
+concept HasEsEQvecFT0M = requires(T collision) {
+  collision.eseQvecFT0MRe();
+  collision.eseQvecFT0MIm();
+};
+
+template <typename T>
+concept HasEsEQvecFV0A = requires(T collision) {
+  collision.eseQvecFV0ARe();
+  collision.eseQvecFV0AIm();
+};
+
+template <typename T>
+concept HasEsEQvecTPCpos = requires(T collision) {
+  collision.eseQvecTPCposRe();
+  collision.eseQvecTPCposIm();
+};
+
+template <typename T>
+concept HasEsEQvecTPCneg = requires(T collision) {
+  collision.eseQvecTPCnegRe();
+  collision.eseQvecTPCnegIm();
+};
+
+template <typename T>
+concept HasEsEQvecTPCtot = requires(T collision) {
+  collision.eseQvecTPCallRe();
+  collision.eseQvecTPCallIm();
+};
+
 /// Get the Ese Q vector choosing your favourite estimator
 /// \param collision is the collision with the Q vector information
 /// \param qvecEst is the chosen Q-vector estimator
@@ -210,37 +253,37 @@ std::array<float, 3> getEseQvec(TCollision const& collision, const int qvecEst)
 {
   switch (qvecEst) {
     case QvecEstimator::FV0A:
-      if constexpr (HasQvecFV0A<TCollision>) {
+      if constexpr (HasEsEQvecFV0A<TCollision>) {
         return std::array<float, 3>{collision.eseQvecFV0ARe(), collision.eseQvecFV0AIm(), collision.sumAmplFV0A()};
       }
       break;
     case QvecEstimator::FT0A:
-      if constexpr (HasQvecFT0A<TCollision>) {
+      if constexpr (HasEsEQvecFT0A<TCollision>) {
         return std::array<float, 3>{collision.eseQvecFT0ARe(), collision.eseQvecFT0AIm(), collision.sumAmplFT0A()};
       }
       break;
     case QvecEstimator::FT0C:
-      if constexpr (HasQvecFT0C<TCollision>) {
+      if constexpr (HasEsEQvecFT0C<TCollision>) {
         return std::array<float, 3>{collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm(), collision.sumAmplFT0C()};
       }
       break;
     case QvecEstimator::FT0M:
-      if constexpr (HasQvecFT0M<TCollision>) {
+      if constexpr (HasEsEQvecFT0M<TCollision>) {
         return std::array<float, 3>{collision.eseQvecFT0MRe(), collision.eseQvecFT0MIm(), collision.sumAmplFT0M()};
       }
       break;
     case QvecEstimator::TPCPos:
-      if constexpr (HasQvecTPCpos<TCollision>) {
+      if constexpr (HasEsEQvecTPCpos<TCollision>) {
         return std::array<float, 3>{collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm(), static_cast<float>(collision.nTrkTPCpos())};
       }
       break;
     case QvecEstimator::TPCNeg:
-      if constexpr (HasQvecTPCneg<TCollision>) {
+      if constexpr (HasEsEQvecTPCneg<TCollision>) {
         return std::array<float, 3>{collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm(), static_cast<float>(collision.nTrkTPCneg())};
       }
       break;
     case QvecEstimator::TPCTot:
-      if constexpr (HasQvecTPCtot<TCollision>) {
+      if constexpr (HasEsEQvecTPCtot<TCollision>) {
         return std::array<float, 3>{collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm(), static_cast<float>(collision.nTrkTPCall())};
       }
       break;