diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx
index 31704c12515..8231ec6858b 100644
--- a/PWGDQ/Core/VarManager.cxx
+++ b/PWGDQ/Core/VarManager.cxx
@@ -635,10 +635,14 @@ void VarManager::SetDefaultVarNames()
   fgVariableUnits[kS11A] = "";
   fgVariableNames[kS12A] = "S_{12}^{A} ";
   fgVariableUnits[kS12A] = "";
+  fgVariableNames[kS21A] = "S_{21}^{A} ";
+  fgVariableUnits[kS21A] = "";
   fgVariableNames[kS13A] = "S_{13}^{A} ";
   fgVariableUnits[kS13A] = "";
   fgVariableNames[kS31A] = "S_{31}^{A} ";
   fgVariableUnits[kS31A] = "";
+  fgVariableNames[kS22A] = "S_{22}^{A} ";
+  fgVariableUnits[kS22A] = "";
   fgVariableNames[kM11REF] = "M_{11}^{REF} ";
   fgVariableUnits[kM11REF] = "";
   fgVariableNames[kM01POI] = "M_{01}^{POI} ";
diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h
index 6def4cb1210..4ac3278b55d 100644
--- a/PWGDQ/Core/VarManager.h
+++ b/PWGDQ/Core/VarManager.h
@@ -518,8 +518,12 @@ class VarManager : public TObject
     kQ23YA,
     kS11A,
     kS12A,
+    kS21A,
     kS13A,
     kS31A,
+    kS22A,
+    kS41A,
+    kS14A,
     kM11REF,
     kM01POI,
     kM1111REF,
@@ -1353,6 +1357,7 @@ void VarManager::FillEvent(T const& event, float* values)
       values[kQ23YA] = event.q23ya();
       values[kS11A] = event.s11a();
       values[kS12A] = event.s12a();
+      values[kS21A] = event.s21a();
       values[kS13A] = event.s13a();
       values[kS31A] = event.s31a();
     }
@@ -3117,8 +3122,12 @@ void VarManager::FillQVectorFromGFW(C const& /*collision*/, A const& compA11, A
   values[kQ23YA] = compA23.imag(); // Only being used by cumulants, no need for normalization
   values[kS11A] = S11A;
   values[kS12A] = S12A;
+  values[kS21A] = S21A;
   values[kS13A] = S13A;
   values[kS31A] = S31A;
+  values[kS22A] = S22A;
+  values[kS14A] = S14A;
+  values[kS41A] = S41A;
 
   // Fill event multiplicities
   values[kMultA] = S10A;
@@ -3303,6 +3312,10 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values)
     values[kM0111POI] = values[kMultDimuons] * (values[kS31A] - 3. * values[kS11A] * values[kS12A] + 2. * values[kS13A]);
     values[kCORR2POI] = (P2 * conj(Q21)).real() / values[kM01POI];
     values[kCORR4POI] = (P2 * Q21 * conj(Q21) * conj(Q21) - P2 * Q21 * conj(Q42) - 2. * values[kS12A] * P2 * conj(Q21) + 2. * P2 * conj(Q23)).real() / values[kM0111POI];
+    values[kM11REF] = values[kS21A] - values[kS12A];
+    values[kM1111REF] = values[kS41A] - 6. * values[kS12A] * values[kS21A] + 8. * values[kS13A] * values[kS11A] + 3. * values[kS22A] - 6. * values[kS14A];
+    values[kCORR2REF] = (norm(Q21) - values[kS12A]) / values[kM11REF];
+    values[kCORR4REF] = (pow(norm(Q21), 2) + norm(Q42) - 2. * (Q42 * conj(Q21) * conj(Q21)).real() + 8. * (Q23 * conj(Q21)).real() - 4. * values[kS12A] * norm(Q21) - 6. * values[kS14A] - 2. * values[kS22A]) / values[kM1111REF];
   }
 }
 
diff --git a/PWGDQ/DataModel/ReducedInfoTables.h b/PWGDQ/DataModel/ReducedInfoTables.h
index 1ef7a3de9e5..cf315e12dd7 100644
--- a/PWGDQ/DataModel/ReducedInfoTables.h
+++ b/PWGDQ/DataModel/ReducedInfoTables.h
@@ -76,6 +76,10 @@ DECLARE_SOA_COLUMN(S11A, s11a, float);         //! Weighted multiplicity (p = 1,
 DECLARE_SOA_COLUMN(S12A, s12a, float);         //! Weighted multiplicity (p = 1, k = 2)
 DECLARE_SOA_COLUMN(S13A, s13a, float);         //! Weighted multiplicity (p = 1, k = 3)
 DECLARE_SOA_COLUMN(S31A, s31a, float);         //! Weighted multiplicity (p = 3, k = 1)
+DECLARE_SOA_COLUMN(S21A, s21a, float);         //! Weighted multiplicity (p = 2, k = 1)
+DECLARE_SOA_COLUMN(S22A, s22a, float);         //! Weighted multiplicity (p = 2, k = 2)
+DECLARE_SOA_COLUMN(S41A, s41a, float);         //! Weighted multiplicity (p = 4, k = 1)
+DECLARE_SOA_COLUMN(S14A, s14a, float);         //! Weighted multiplicity (p = 1, k = 4)
 DECLARE_SOA_COLUMN(CORR2REF, corr2ref, float); //!  Ref Flow correlator <2>
 DECLARE_SOA_COLUMN(CORR4REF, corr4ref, float); //!  Ref Flow correlator <4>
 DECLARE_SOA_COLUMN(M11REF, m11ref, float);     //!  Weighted multiplicity of <<2>> for reference flow
@@ -110,7 +114,7 @@ DECLARE_SOA_TABLE(ReducedEventsQvector, "AOD", "REQVECTOR", //!    Event Q-vecto
 
 DECLARE_SOA_TABLE(ReducedEventsQvectorExtra, "AOD", "REQVECTOREXTRA", //!    Event Q-vector extra information
                   reducedevent::Q42XA, reducedevent::Q42YA, reducedevent::Q23XA, reducedevent::Q23YA,
-                  reducedevent::S11A, reducedevent::S12A, reducedevent::S13A, reducedevent::S31A);
+                  reducedevent::S11A, reducedevent::S12A, reducedevent::S21A, reducedevent::S13A, reducedevent::S31A);
 
 DECLARE_SOA_TABLE(ReducedEventsQvectorCentr, "AOD", "REQVECTORCTR", //!    Event Q-vector information from central framework
                   qvec::QvecFT0ARe, qvec::QvecFT0AIm, qvec::QvecFT0CRe, qvec::QvecFT0CIm, qvec::QvecFT0MRe, qvec::QvecFT0MIm, qvec::QvecFV0ARe, qvec::QvecFV0AIm, qvec::QvecBPosRe, qvec::QvecBPosIm, qvec::QvecBNegRe, qvec::QvecBNegIm,
@@ -530,6 +534,11 @@ DECLARE_SOA_COLUMN(R2SP, r2sp, float);                                   //! Eve
 DECLARE_SOA_COLUMN(R2EP, r2ep, float);                                   //! Event plane resolution for EP method
 DECLARE_SOA_COLUMN(CORR2POI, corr2poi, float);                           //! POI FLOW CORRELATOR <2'>
 DECLARE_SOA_COLUMN(CORR4POI, corr4poi, float);                           //! POI FLOW CORRELATOR <4'>
+DECLARE_SOA_COLUMN(CORR2REF, corr2ref, float);                           //! POI FLOW CORRELATOR <2> (by dimuons)
+DECLARE_SOA_COLUMN(CORR4REF, corr4ref, float);                           //! POI FLOW CORRELATOR <4> (by dimuons)
+DECLARE_SOA_COLUMN(M11REF, m11ref, float);                               //!  Weighted multiplicity of <<2>> for reference flow (by dimuons)
+DECLARE_SOA_COLUMN(M1111REF, m1111ref, float);                           //!  Weighted multiplicity of <<4>> for reference flow (by dimuons)
+DECLARE_SOA_COLUMN(MultA, multa, float);                                 //!  Multiplicity A of reference flow (by dimuons)
 DECLARE_SOA_COLUMN(M01POI, m01poi, float);                               //! POI event weight for <2'>
 DECLARE_SOA_COLUMN(M0111POI, m0111poi, float);                           //! POI event weight for <4'>
 DECLARE_SOA_COLUMN(MultDimuons, multdimuons, int);                       //! Dimuon multiplicity
@@ -586,6 +595,14 @@ DECLARE_SOA_TABLE(DileptonsFlow, "AOD", "RTDILEPTONFLOW", //!
                   reducedpair::Cos2DeltaPhi,
                   reducedpair::Cos3DeltaPhi);
 
+DECLARE_SOA_TABLE(RefFlowDimuons, "AOD", "RTREFFLOWDIMUON", //!
+                  reducedpair::CORR2REF,
+                  reducedpair::CORR4REF,
+                  reducedpair::M11REF,
+                  reducedpair::M1111REF,
+                  reducedpair::CentFT0C,
+                  reducedpair::MultA);
+
 // Dilepton collision information (joined with DileptonsExtra) allowing to connect different tables (cross PWGs)
 DECLARE_SOA_TABLE(DileptonsInfo, "AOD", "RTDILEPTONINFO",
                   reducedpair::CollisionId, collision::PosX, collision::PosY, collision::PosZ);
@@ -631,6 +648,7 @@ using Dimuon = Dimuons::iterator;
 using DielectronExtra = DielectronsExtra::iterator;
 using DimuonExtra = DimuonsExtra::iterator;
 using DileptonFlow = DileptonsFlow::iterator;
+using RefFlowDimuon = RefFlowDimuons::iterator;
 using DileptonInfo = DileptonsInfo::iterator;
 using DimuonAll = DimuonsAll::iterator;
 
diff --git a/PWGDQ/Tasks/dqFlow.cxx b/PWGDQ/Tasks/dqFlow.cxx
index 3eb3cde054a..196eeacb665 100644
--- a/PWGDQ/Tasks/dqFlow.cxx
+++ b/PWGDQ/Tasks/dqFlow.cxx
@@ -451,7 +451,7 @@ struct DQEventQvector {
     // Fill the tree for the reduced event table with Q vector quantities
     if (fEventCut->IsSelected(VarManager::fgValues)) {
       eventQvector(VarManager::fgValues[VarManager::kQ1X0A], VarManager::fgValues[VarManager::kQ1Y0A], VarManager::fgValues[VarManager::kQ1X0B], VarManager::fgValues[VarManager::kQ1Y0B], VarManager::fgValues[VarManager::kQ1X0C], VarManager::fgValues[VarManager::kQ1Y0C], VarManager::fgValues[VarManager::kQ2X0A], VarManager::fgValues[VarManager::kQ2Y0A], VarManager::fgValues[VarManager::kQ2X0B], VarManager::fgValues[VarManager::kQ2Y0B], VarManager::fgValues[VarManager::kQ2X0C], VarManager::fgValues[VarManager::kQ2Y0C], VarManager::fgValues[VarManager::kMultA], VarManager::fgValues[VarManager::kMultB], VarManager::fgValues[VarManager::kMultC], VarManager::fgValues[VarManager::kQ3X0A], VarManager::fgValues[VarManager::kQ3Y0A], VarManager::fgValues[VarManager::kQ3X0B], VarManager::fgValues[VarManager::kQ3Y0B], VarManager::fgValues[VarManager::kQ3X0C], VarManager::fgValues[VarManager::kQ3Y0C], VarManager::fgValues[VarManager::kQ4X0A], VarManager::fgValues[VarManager::kQ4Y0A], VarManager::fgValues[VarManager::kQ4X0B], VarManager::fgValues[VarManager::kQ4Y0B], VarManager::fgValues[VarManager::kQ4X0C], VarManager::fgValues[VarManager::kQ4Y0C]);
-      eventQvectorExtra(VarManager::fgValues[VarManager::kQ42XA], VarManager::fgValues[VarManager::kQ42YA], VarManager::fgValues[VarManager::kQ23XA], VarManager::fgValues[VarManager::kQ23YA], VarManager::fgValues[VarManager::kS11A], VarManager::fgValues[VarManager::kS12A], VarManager::fgValues[VarManager::kS13A], VarManager::fgValues[VarManager::kS31A]);
+      eventQvectorExtra(VarManager::fgValues[VarManager::kQ42XA], VarManager::fgValues[VarManager::kQ42YA], VarManager::fgValues[VarManager::kQ23XA], VarManager::fgValues[VarManager::kQ23YA], VarManager::fgValues[VarManager::kS11A], VarManager::fgValues[VarManager::kS12A], VarManager::fgValues[VarManager::kS21A], VarManager::fgValues[VarManager::kS13A], VarManager::fgValues[VarManager::kS31A]);
       eventRefFlow(VarManager::fgValues[VarManager::kM11REF], VarManager::fgValues[VarManager::kM1111REF], VarManager::fgValues[VarManager::kCORR2REF], VarManager::fgValues[VarManager::kCORR4REF], centrality);
     }
 
diff --git a/PWGDQ/Tasks/tableReader.cxx b/PWGDQ/Tasks/tableReader.cxx
index 3838f40380c..2a09780563f 100644
--- a/PWGDQ/Tasks/tableReader.cxx
+++ b/PWGDQ/Tasks/tableReader.cxx
@@ -756,6 +756,7 @@ struct AnalysisSameEventPairing {
   Produces<aod::DimuonsExtra> dimuonExtraList;
   Produces<aod::DimuonsAll> dimuonAllList;
   Produces<aod::DileptonFlow> dileptonFlowList;
+  Produces<aod::RefFlowDimuon> refFlowDimuonList;
   Produces<aod::DileptonsInfo> dileptonInfoList;
   float mMagField = 0.0;
   o2::parameters::GRPMagField* grpmag = nullptr;
@@ -1108,6 +1109,7 @@ struct AnalysisSameEventPairing {
 
       if constexpr (eventHasQvector) {
         dileptonFlowList(VarManager::fgValues[VarManager::kU2Q2], VarManager::fgValues[VarManager::kU3Q3], VarManager::fgValues[VarManager::kCos2DeltaPhi], VarManager::fgValues[VarManager::kCos3DeltaPhi]);
+        refFlowDimuonList(VarManager::fgValues[VarManager::kCORR2REF], VarManager::fgValues[VarManager::kCORR4REF], VarManager::fgValues[VarManager::kM11REF], VarManager::fgValues[VarManager::kM1111REF], VarManager::fgValues[VarManager::kCentFT0C], VarManager::fgValues[VarManager::kMultA]);
       }
 
       int iCut = 0;
@@ -1253,14 +1255,14 @@ struct AnalysisSameEventPairing {
     VarManager::FillEvent<gkEventFillMapWithCovQvector>(event, VarManager::fgValues);
     runSameEventPairing<true, VarManager::kDecayToEE, gkEventFillMapWithCovQvector, gkTrackFillMap>(event, tracks, tracks);
   }
-  void processVnDecayToMuMuSkimmed(soa::Filtered<MyEventsVtxCovSelectedQvector>::iterator const& event, soa::Filtered<MyMuonTracksSelected> const& muons)
+  void processVnDecayToMuMuSkimmed(soa::Filtered<MyEventsVtxCovSelectedQvectorExtra>::iterator const& event, soa::Filtered<MyMuonTracksSelectedWithColl> const& muons)
   {
     // Reset the fValues array
     VarManager::ResetValues(0, VarManager::kNVars);
-    VarManager::FillEvent<gkEventFillMapWithCovQvector>(event, VarManager::fgValues);
-    runSameEventPairing<true, VarManager::kDecayToMuMu, gkEventFillMapWithCovQvector, gkMuonFillMap>(event, muons, muons);
+    VarManager::FillEvent<gkEventFillMapWithCovQvectorExtra>(event, VarManager::fgValues);
+    runSameEventPairing<true, VarManager::kDecayToMuMu, gkEventFillMapWithCovQvectorExtra, gkMuonFillMap>(event, muons, muons);
   }
-  void processVnDecayToMuMuSkimmedWithWeights(soa::Filtered<MyEventsVtxCovSelectedQvectorExtra>::iterator const& event, soa::Filtered<MyMuonTracksSelected> const& muons)
+  void processVnDecayToMuMuSkimmedWithWeights(soa::Filtered<MyEventsVtxCovSelectedQvectorExtra>::iterator const& event, soa::Filtered<MyMuonTracksSelectedWithColl> const& muons)
   {
     // Reset the fValues array
     VarManager::ResetValues(0, VarManager::kNVars);
diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx
index f34b3d73a51..ae4905b2636 100644
--- a/PWGDQ/Tasks/tableReader_withAssoc.cxx
+++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx
@@ -790,6 +790,7 @@ struct AnalysisSameEventPairing {
   Produces<aod::DimuonsExtra> dimuonsExtraList;
   Produces<aod::DimuonsAll> dimuonAllList;
   Produces<aod::DileptonFlow> dileptonFlowList;
+  Produces<aod::RefFlowDimuon> refFlowDimuonList;
   Produces<aod::DileptonsInfo> dileptonInfoList;
 
   o2::base::MatLayerCylSet* fLUT = nullptr;
@@ -1250,6 +1251,7 @@ struct AnalysisSameEventPairing {
 
         if constexpr (eventHasQvector) {
           dileptonFlowList(VarManager::fgValues[VarManager::kU2Q2], VarManager::fgValues[VarManager::kU3Q3], VarManager::fgValues[VarManager::kCos2DeltaPhi], VarManager::fgValues[VarManager::kCos3DeltaPhi]);
+          refFlowDimuonList(VarManager::fgValues[VarManager::kCORR2REF], VarManager::fgValues[VarManager::kCORR4REF], VarManager::fgValues[VarManager::kM11REF], VarManager::fgValues[VarManager::kM1111REF], VarManager::fgValues[VarManager::kCentFT0C], VarManager::fgValues[VarManager::kMultA]);
         }
 
         // Fill histograms