diff --git a/PWGCF/Flow/Tasks/flowEventPlane.cxx b/PWGCF/Flow/Tasks/flowEventPlane.cxx
index 502da4810b2..655574c9bde 100644
--- a/PWGCF/Flow/Tasks/flowEventPlane.cxx
+++ b/PWGCF/Flow/Tasks/flowEventPlane.cxx
@@ -16,6 +16,7 @@
 #include "PWGLF/DataModel/LFStrangenessTables.h"
 
 #include "Common/CCDB/EventSelectionParams.h"
+#include "Common/Core/RecoDecay.h"
 #include "Common/DataModel/Centrality.h"
 #include "Common/DataModel/CollisionAssociationTables.h"
 #include "Common/DataModel/EventSelection.h"
@@ -24,8 +25,11 @@
 #include "Common/DataModel/PIDResponseTPC.h"
 #include "Common/DataModel/TrackSelectionTables.h"
 
+#include <CCDB/BasicCCDBManager.h>
+#include <CommonConstants/MathConstants.h>
 #include <CommonConstants/PhysicsConstants.h>
 #include <Framework/ASoA.h>
+#include <Framework/ASoAHelpers.h>
 #include <Framework/AnalysisDataModel.h>
 #include <Framework/AnalysisHelpers.h>
 #include <Framework/AnalysisTask.h>
@@ -34,12 +38,19 @@
 #include <Framework/HistogramSpec.h>
 #include <Framework/InitContext.h>
 #include <Framework/OutputObjHeader.h>
+#include <Framework/SliceCache.h>
 #include <Framework/runDataProcessing.h>
 
+#include <TH2.h>
+#include <THnSparse.h>
+#include <TList.h>
+#include <TProfile.h>
+
 #include <array>
 #include <cmath>
-#include <cstdint>
+#include <map>
 #include <string>
+#include <string_view>
 #include <vector>
 
 using namespace o2;
@@ -50,7 +61,7 @@ using namespace o2::constants::math;
 
 namespace o2::aod
 {
-namespace colsp
+namespace colspcalib
 {
 DECLARE_SOA_COLUMN(RunNumber, runNumber, int);
 DECLARE_SOA_COLUMN(Timestamp, timestamp, uint64_t);
@@ -62,75 +73,85 @@ DECLARE_SOA_COLUMN(ZnaEnergyCommon, znaEnergyCommon, float);
 DECLARE_SOA_COLUMN(ZncEnergyCommon, zncEnergyCommon, float);
 DECLARE_SOA_COLUMN(ZnaEnergy, znaEnergy, float[4]);
 DECLARE_SOA_COLUMN(ZncEnergy, zncEnergy, float[4]);
-} // namespace colsp
-DECLARE_SOA_TABLE(ColSps, "AOD", "COLSPS", o2::soa::Index<>,
-                  colsp::RunNumber,
-                  colsp::Timestamp,
-                  colsp::Cent,
-                  colsp::Vx,
-                  colsp::Vy,
-                  colsp::Vz,
-                  colsp::ZnaEnergyCommon,
-                  colsp::ZncEnergyCommon,
-                  colsp::ZnaEnergy,
-                  colsp::ZncEnergy);
-
-using ColSp = ColSps::iterator;
+} // namespace colspcalib
+DECLARE_SOA_TABLE(ColSpCalib, "AOD", "COLSPCALIB", o2::soa::Index<>,
+                  colspcalib::RunNumber,
+                  colspcalib::Timestamp,
+                  colspcalib::Cent,
+                  colspcalib::Vx,
+                  colspcalib::Vy,
+                  colspcalib::Vz,
+                  colspcalib::ZnaEnergyCommon,
+                  colspcalib::ZncEnergyCommon,
+                  colspcalib::ZnaEnergy,
+                  colspcalib::ZncEnergy);
+
+namespace colspext
+{
+DECLARE_SOA_COLUMN(SelColFlag, selColFlag, bool);
+DECLARE_SOA_COLUMN(Xa, xa, float);
+DECLARE_SOA_COLUMN(Ya, ya, float);
+DECLARE_SOA_COLUMN(Xc, xc, float);
+DECLARE_SOA_COLUMN(Yc, yc, float);
+} // namespace colspext
+DECLARE_SOA_TABLE(ColSPExt, "AOD", "COLSPSEXT", o2::soa::Index<>,
+                  colspext::SelColFlag,
+                  colspext::Xa,
+                  colspext::Ya,
+                  colspext::Xc,
+                  colspext::Yc);
 
 namespace tracksid
 {
-DECLARE_SOA_INDEX_COLUMN(ColSp, colSp);
-DECLARE_SOA_COLUMN(Px, px, float);
-DECLARE_SOA_COLUMN(Py, py, float);
-DECLARE_SOA_COLUMN(Pz, pz, float);
-DECLARE_SOA_COLUMN(Charge, charge, int8_t);
-DECLARE_SOA_COLUMN(TpcInfo, tpcInfo, float[4]);
-DECLARE_SOA_COLUMN(TofInfo, tofInfo, float[4]);
-DECLARE_SOA_COLUMN(Dca, dca, float[2]);
+DECLARE_SOA_COLUMN(IsCharged, isCharged, bool);
+DECLARE_SOA_COLUMN(IsPion, isPion, bool);
+DECLARE_SOA_COLUMN(IsKaon, isKaon, bool);
+DECLARE_SOA_COLUMN(IsProton, isProton, bool);
 } // namespace tracksid
 DECLARE_SOA_TABLE(TracksId, "AOD", "TRACKSID", o2::soa::Index<>,
-                  tracksid::ColSpId,
-                  tracksid::Px,
-                  tracksid::Py,
-                  tracksid::Pz,
-                  tracksid::Charge,
-                  tracksid::TpcInfo,
-                  tracksid::TofInfo,
-                  tracksid::Dca);
-using TrackId = TracksId::iterator;
-
-namespace v0track
-{
-DECLARE_SOA_INDEX_COLUMN(ColSp, colSp);
-DECLARE_SOA_COLUMN(Px, px, float);
-DECLARE_SOA_COLUMN(Py, py, float);
-DECLARE_SOA_COLUMN(Pz, pz, float);
-DECLARE_SOA_COLUMN(CosPA, cosPA, float);
-DECLARE_SOA_COLUMN(Ctau, ctau, float);
-DECLARE_SOA_COLUMN(Mass, mass, float);
-DECLARE_SOA_COLUMN(Species, species, uint8_t);
-DECLARE_SOA_COLUMN(DcaDau, dcaDau, float[2]);
-DECLARE_SOA_COLUMN(TpcDau, tpcDau, float[2]);
-} // namespace v0track
-DECLARE_SOA_TABLE(V0Tracks, "AOD", "V0TRACKS", o2::soa::Index<>,
-                  v0track::ColSpId,
-                  v0track::Px,
-                  v0track::Py,
-                  v0track::Pz,
-                  v0track::CosPA,
-                  v0track::Ctau,
-                  v0track::Mass,
-                  v0track::Species,
-                  v0track::DcaDau,
-                  v0track::TpcDau);
-using V0Track = V0Tracks::iterator;
+                  tracksid::IsCharged,
+                  tracksid::IsPion,
+                  tracksid::IsKaon,
+                  tracksid::IsProton);
 } // namespace o2::aod
 
-enum TrackType {
-  kCharged = 0,
-  kPi,
+enum GainClibCorr {
+  kGainCalibA = 0,
+  kGainCalibC,
+  kNGainCalib
+};
+
+enum CorrectionType {
+  kFineCorr = 0,
+  kCoarseCorr,
+  kNCorr
+};
+
+enum CollisionParameterType {
+  kCent = 0,
+  kVx,
+  kVy,
+  kVz
+};
+
+enum ZDCXYType {
+  kXa = 0,
+  kYa,
+  kXc,
+  kYc,
+  kXYAC
+};
+
+enum ParticleType {
+  kPi = 0,
   kKa,
-  kPr
+  kPr,
+  kNPart
+};
+
+enum ResoType {
+  kPhi0 = 0,
+  kKStar
 };
 
 enum V0Type {
@@ -139,11 +160,134 @@ enum V0Type {
   kAntiLambda
 };
 
-struct FlowEventPlane {
+struct SpCalibTableProducer {
+  // Table producer
+  Produces<aod::ColSpCalib> colSpCalibTable;
+
+  // Configurables
+  // Collisions
+  Configurable<float> cMinZVtx{"cMinZVtx", -10.0, "Min VtxZ cut"};
+  Configurable<float> cMaxZVtx{"cMaxZVtx", 10.0, "Max VtxZ cut"};
+  Configurable<float> cMinCent{"cMinCent", 0., "Minumum Centrality"};
+  Configurable<float> cMaxCent{"cMaxCent", 100.0, "Maximum Centrality"};
+  Configurable<bool> cSel8Trig{"cSel8Trig", true, "Sel8 (T0A + T0C) Selection Run3"};
+  Configurable<bool> cPileupReject{"cPileupReject", true, "Pileup rejection"};
+  Configurable<bool> cZVtxTimeDiff{"cZVtxTimeDiff", true, "z-vtx time diff selection"};
+  Configurable<bool> cIsGoodITSLayers{"cIsGoodITSLayers", true, "Good ITS Layers All"};
+  Configurable<float> cMinOccupancy{"cMinOccupancy", 0, "Minimum FT0C Occupancy"};
+  Configurable<float> cMaxOccupancy{"cMaxOccupancy", 1e6, "Maximum FT0C Occupancy"};
+
+  // Histogram Registry.
+  HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
+
+  // Run number
+  uint64_t runNum = 0, timestamp = 0;
+  float mult = 0., cent = 0., posX = 0., posY = 0., posZ = 0.;
+
+  void init(InitContext const&)
+  {
+    // Histogram collision
+    histos.add("hCent", "Centrality", kTH1F, {{100, 0., 100., "FT0C%"}});
+    histos.add("hVz", "V_{z}", kTH1F, {{100, -10., 10., "V_{z}(cm)"}});
+  }
+
+  template <typename C>
+  bool selCollision(C const& col)
+  {
+    if (col.posZ() <= cMinZVtx || col.posZ() >= cMaxZVtx) { // VtxZ selection
+      return false;
+    }
+
+    if (cSel8Trig && !col.sel8()) { // Sel8 selection
+      return false;
+    }
+
+    cent = col.centFT0C();
+    if (cent <= cMinCent || cent >= cMaxCent) { // Centrality selection
+      return false;
+    }
+
+    if (col.ft0cOccupancyInTimeRange() < cMinOccupancy || col.ft0cOccupancyInTimeRange() > cMaxOccupancy) { // Occupancy cut
+      return false;
+    }
+
+    if (cPileupReject && !col.selection_bit(aod::evsel::kNoSameBunchPileup)) { // Pile-up rejection
+      return false;
+    }
+
+    if (cZVtxTimeDiff && !col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { // ZvtxFT0 vs PV
+      return false;
+    }
+
+    if (cIsGoodITSLayers && !col.selection_bit(aod::evsel::kIsGoodITSLayersAll)) { // All ITS layer active
+      return false;
+    }
+
+    // Set Multiplicity
+    mult = col.multTPC();
+
+    return true;
+  }
+
+  template <typename B, typename C>
+  void analyzeCollision(B const& bc, C const& collision)
+  {
+    // Check zdc
+    if (!bc.has_zdc()) {
+      return;
+    }
+
+    // Event selection
+    if (!selCollision(collision)) {
+      return;
+    }
+    posX = collision.posX();
+    posY = collision.posY();
+    posZ = collision.posZ();
+
+    auto zdc = bc.zdc();
+    auto znaEnergy = zdc.energySectorZNA();
+    auto zncEnergy = zdc.energySectorZNC();
+    auto znaEnergyCommon = zdc.energyCommonZNA();
+    auto zncEnergyCommon = zdc.energyCommonZNC();
+
+    // check energy deposits
+    if (znaEnergyCommon <= 0 || zncEnergyCommon <= 0 || znaEnergy[0] <= 0 || znaEnergy[1] <= 0 || znaEnergy[2] <= 0 || znaEnergy[3] <= 0 || zncEnergy[0] <= 0 || zncEnergy[1] <= 0 || zncEnergy[2] <= 0 || zncEnergy[3] <= 0) {
+      return;
+    }
+
+    // Fill collision table
+    histos.fill(HIST("hCent"), cent);
+    histos.fill(HIST("hVz"), posZ);
+    colSpCalibTable(runNum, timestamp, cent, posX, posY, posZ, znaEnergyCommon, zncEnergyCommon, znaEnergy.data(), zncEnergy.data());
+
+    // Done
+    return;
+  }
+
+  using BCsRun3 = soa::Join<aod::BCsWithTimestamps, aod::Run3MatchedToBCSparse>;
+  using CollisionsRun3 = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::CentFT0Ms, aod::CentFV0As, aod::TPCMults, aod::PVMults, aod::MultsGlobal, aod::MultsExtra>;
+
+  void processDummy(CollisionsRun3::iterator const&) {}
+
+  PROCESS_SWITCH(SpCalibTableProducer, processDummy, "Dummy process", true);
+
+  void processFEP(CollisionsRun3::iterator const& collision, BCsRun3 const&, aod::Zdcs const&)
+  {
+    // Get bunch crossing
+    auto bc = collision.template foundBC_as<BCsRun3>();
+    runNum = collision.template foundBC_as<BCsRun3>().runNumber();
+    timestamp = collision.template foundBC_as<BCsRun3>().timestamp();
+    analyzeCollision(bc, collision);
+  }
+
+  PROCESS_SWITCH(SpCalibTableProducer, processFEP, "Flow Event Plane Table Producer", false);
+};
+
+struct SpectatorPlaneTableProducer {
   // Table producer
-  Produces<aod::ColSp> colSpTable;
+  Produces<aod::ColSPExt> colSPExtTable;
   Produces<aod::TracksId> tracksIdTable;
-  Produces<aod::V0Tracks> v0sTable;
 
   // Configurables
   // Collisions
@@ -158,9 +302,35 @@ struct FlowEventPlane {
   Configurable<float> cMinOccupancy{"cMinOccupancy", 0, "Minimum FT0C Occupancy"};
   Configurable<float> cMaxOccupancy{"cMaxOccupancy", 1e6, "Maximum FT0C Occupancy"};
 
+  // Coarse binning factor
+  Configurable<int> cAxisCBF{"cAxisCBF", 10, "Coarse Bin Factor"};
+
+  // Cent Vx Vy Vz Bins
+  Configurable<int> cAxisCentBins{"cAxisCentBins", 50, "NBins Centrality"};
+  Configurable<int> cAxisVxyBins{"cAxisVxyBins", 50, "NBins Vx Vy"};
+  Configurable<int> cAxisVzBins{"cAxisVzBins", 50, "NBins Vz"};
+  Configurable<float> cAxisVxMin{"cAxisVxMin", -0.01, "Vx Min"};
+  Configurable<float> cAxisVxMax{"cAxisVxMax", 0.01, "Vx Max"};
+  Configurable<float> cAxisVyMin{"cAxisVyMin", -0.01, "Vy Min"};
+  Configurable<float> cAxisVyMax{"cAxisVyMax", 0.01, "Vy Max"};
+  Configurable<bool> cRecentVxVy{"cRecentVxVy", true, "Vx / Vy recentering"};
+
+  // Corrections
+  Configurable<bool> cLoadCorrection{"cLoadCorrection", true, "Load corrections"};
+  Configurable<bool> cDoGainCalib{"cDoGainCalib", true, "Gain Calib Flag"};
+  Configurable<bool> cUseAlphaZDC{"cUseAlphaZDC", true, "Use Alpha ZDC"};
+  Configurable<bool> cApplyRecentCorr{"cApplyRecentCorr", true, "Apply recentering"};
+  Configurable<std::vector<int>> cCorrFlagVector{"cCorrFlagVector", {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, "Correction Flag"};
+
+  // CCDB
+  Configurable<std::string> cCcdbUrl{"cCcdbUrl", "http://ccdb-test.cern.ch:8080", "url of ccdb"};
+  Configurable<std::string> cCcdbPath{"cCcdbPath", "Users/y/ypatley/DFOO", "Path for ccdb-object"};
+  Configurable<int64_t> nolaterthan{"nolaterthan", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"};
+
   // Tracks
   Configurable<float> cTrackMinPt{"cTrackMinPt", 0.1, "p_{T} minimum"};
   Configurable<float> cTrackMaxPt{"cTrackMaxPt", 10.0, "p_{T} maximum"};
+  Configurable<int> cNEtaBins{"cNEtaBins", 7, "# of eta bins"};
   Configurable<float> cTrackEtaCut{"cTrackEtaCut", 0.8, "Pseudorapidity cut"};
   Configurable<bool> cTrackGlobal{"cTrackGlobal", true, "Global Track"};
   Configurable<float> cTrackDcaXYCut{"cTrackDcaXYCut", 0.1, "DcaXY Cut"};
@@ -169,60 +339,141 @@ struct FlowEventPlane {
   // Track PID
   Configurable<float> cTpcElRejCutMin{"cTpcElRejCutMin", -3., "Electron Rejection Cut Minimum"};
   Configurable<float> cTpcElRejCutMax{"cTpcElRejCutMax", 5., "Electron Rejection Cut Maximum"};
-  Configurable<float> cTpcNSigmaCut{"cTpcNSigmaCut", 5, "TPC NSigma Cut"};
+  Configurable<float> cTpcNSigmaCut{"cTpcNSigmaCut", 2, "TPC NSigma Cut"};
   Configurable<float> cTpcRejCut{"cTpcRejCut", 3, "TPC Rej Cut"};
-  Configurable<float> cTofNSigmaCut{"cTofNSigmaCut", 5, "TOF NSigma Cut"};
+  Configurable<float> cTofNSigmaCut{"cTofNSigmaCut", 2, "TOF NSigma Cut"};
   Configurable<float> cTofRejCut{"cTofRejCut", 3, "TOF Rej Cut"};
+  Configurable<float> cPionPtCut{"cPionPtCut", 0.6, "Pion TPC pT cutoff"};
+  Configurable<float> cKaonPtCut{"cKaonPtCut", 0.6, "Kaon TPC pT cutoff"};
+  Configurable<float> cProtonPtCut{"cProtonPtCut", 1.1, "Proton TPC pT cutoff"};
 
-  // V0s
-  Configurable<int> cV0TypeSelection{"cV0TypeSelection", 1, "V0 Type Selection"};
-  Configurable<float> cMinDcaProtonToPV{"cMinDcaProtonToPV", 0.02, "Minimum Proton DCAr to PV"};
-  Configurable<float> cMinDcaPionToPV{"cMinDcaPionToPV", 0.06, "Minimum Pion DCAr to PV"};
-  Configurable<float> cDcaV0Dau{"cDcaV0Dau", 1., "DCA between V0 daughters"};
-  Configurable<float> cMinV0Radius{"cMinV0Radius", 0.5, "Minimum V0 radius from PV"};
-  Configurable<float> cV0CTau{"cV0CTau", 50.0, "Decay length selection"};
-  Configurable<float> cV0CosPA{"cV0CosPA", 0.95, "CosPA selection"};
-  Configurable<float> cArmPodSel{"cArmPodSel", 0.2, "Armentros-Podolanski Selection for K0S"};
-  Configurable<float> cV0EtaCut{"cV0EtaCut", 0.8, "V0 eta cut"};
+  // Initialize CCDB Service
+  Service<o2::ccdb::BasicCCDBManager> ccdbService;
 
-  // Histogram Registry.
+  // Histogram registry: an object to hold your histograms
   HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
 
+  // Global objects
+  const float zdcDenThrs = 1e-4;
+  float cent = 0., mult = 0.;
+  float posX = 0., posY = 0., posZ = 0.;
+  std::vector<std::vector<std::string>> vCoarseCorrHistNames = {
+    {"hXZNAVsCentVxVyVz"},
+    {"hYZNAVsCentVxVyVz"},
+    {"hXZNCVsCentVxVyVz"},
+    {"hYZNCVsCentVxVyVz"}};
+  std::vector<std::vector<std::string>> vFineCorrHistNames = {
+    {"hXZNAVsCent", "hXZNAVsVx", "hXZNAVsVy", "hXZNAVsVz"},
+    {"hYZNAVsCent", "hYZNAVsVx", "hYZNAVsVy", "hYZNAVsVz"},
+    {"hXZNCVsCent", "hXZNCVsVx", "hXZNCVsVy", "hXZNCVsVz"},
+    {"hYZNCVsCent", "hYZNCVsVx", "hYZNCVsVy", "hYZNCVsVz"}};
+  std::map<CorrectionType, std::vector<std::vector<std::string>>> corrTypeHistNameMap = {{kFineCorr, vFineCorrHistNames}, {kCoarseCorr, vCoarseCorrHistNames}};
+
+  // Container for histograms
+  struct CorrectionHistContainer {
+    TProfile* hVx;
+    TProfile* hVy;
+    std::array<TH2F*, 2> hGainCalib;
+    std::array<std::array<std::array<THnSparseF*, 1>, 4>, 14> vCoarseCorrHist;
+    std::array<std::array<std::array<TProfile*, 4>, 4>, 14> vFineCorrHist;
+  } CorrectionHistContainer;
+
   // Run number
-  uint64_t runNum = 0, timestamp = 0;
-  int64_t colIdx = 0;
-  float mult = 0., cent = 0., posX = 0., posY = 0., posZ = 0.;
+  int cRunNum = 0, lRunNum = 0;
 
   void init(InitContext const&)
   {
-    // Axis collision
-    const AxisSpec axisCent{100, 0., 100, "FT0C%"};
+    // Set CCDB url
+    ccdbService->setURL(cCcdbUrl.value);
+    ccdbService->setCaching(true);
+    ccdbService->setLocalObjectValidityChecking();
+    ccdbService->setCreatedNotAfter(nolaterthan.value);
 
-    const AxisSpec axisTrackPt(40, 0, 4, "p_{T} (GeV/#it{c})");
-    const AxisSpec axisTrackDcaZ(220, -1.1, 1.1, "dca_{Z} (cm)");
-    const AxisSpec axisTrackDcaXY(80, -0.2, 0.2, "dca_{XY} (cm)");
-    const AxisSpec axisTrackdEdx(360, 20, 200, "#frac{dE}{dx}");
-    const AxisSpec axisTrackTofSignal(240, 0, 1.2, "#beta");
+    // Define axes
+    const AxisSpec axisZDCEnergy{500, 0, 500, "ZD[AC] Signal"};
 
-    const AxisSpec axisAlpha(40, -1, 1, "#alpha");
-    const AxisSpec axisQtarm(40, 0, 0.4, "q_{T}");
-
-    // Histogram collision
-    histos.add("hCent", "Centrality", kTH1F, {axisCent});
-
-    // Tracks QA
-    histos.add("QA/hDcaXY", "DCA_{XY} vs pT", kTH2F, {axisTrackPt, axisTrackDcaXY});
-    histos.add("QA/hDcaZ", "DCA_{Z} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaZ});
-    histos.add("QA/hdEdX", "dE/dx vs pT", kTH2F, {axisTrackPt, axisTrackdEdx});
-    histos.add("QA/hTOFSignal", "#beta_{TOF} vs p_{T}", kTH2F, {axisTrackPt, axisTrackTofSignal});
-
-    // V0s QA
-    histos.add("QA/h2f_armpod_before_sel", "Armentros-Podolanski Plot", kTH2F, {axisAlpha, axisQtarm});
-    histos.add("QA/h2f_armpod_K0s", "Armentros-Podolanski Plot", kTH2F, {axisAlpha, axisQtarm});
-    histos.add("QA/h2f_armpod_Lambda", "Armentros-Podolanski Plot", kTH2F, {axisAlpha, axisQtarm});
-    histos.add("QA/h2f_armpod_AntiLambda", "Armentros-Podolanski Plot", kTH2F, {axisAlpha, axisQtarm});
+    const AxisSpec axisCent{100, 0., 100, "FT0C%"};
+    const AxisSpec axisVx{cAxisVxyBins, cAxisVxMin, cAxisVxMax, "V_{X}(cm)"};
+    const AxisSpec axisVy{cAxisVxyBins, cAxisVyMin, cAxisVyMax, "V_{Y}(cm)"};
+    const AxisSpec axisVz{cAxisVzBins, cMinZVtx, cMaxZVtx, "V_{Z}(cm)"};
+
+    const AxisSpec axisCoarseCent{cAxisCentBins / cAxisCBF, cMinCent, cMaxCent, "FT0C%"};
+    const AxisSpec axisCoarseVx{cAxisVxyBins / cAxisCBF, cAxisVxMin, cAxisVxMax, "V_{x}"};
+    const AxisSpec axisCoarseVy{cAxisVxyBins / cAxisCBF, cAxisVyMin, cAxisVyMax, "V_{y}"};
+    const AxisSpec axisCoarseVz{cAxisVzBins / cAxisCBF, cMinZVtx, cMaxZVtx, "V_{z}"};
+
+    const AxisSpec axisFineCent{cAxisCentBins, cMinCent, cMaxCent, "FT0C%"};
+    const AxisSpec axisFineVx{cAxisVxyBins, cAxisVxMin, cAxisVxMax, "V_{x}"};
+    const AxisSpec axisFineVy{cAxisVxyBins, cAxisVyMin, cAxisVyMax, "V_{x}"};
+    const AxisSpec axisFineVz{cAxisVzBins, cMinZVtx, cMaxZVtx, "V_{z}"};
+
+    const AxisSpec axisXa{40, -1, 1, "X^{ZNA}_{1}"};
+    const AxisSpec axisYa{40, -1, 1, "Y^{ZNA}_{1}"};
+    const AxisSpec axisXc{40, -1, 1, "X^{ZNC}_{1}"};
+    const AxisSpec axisYc{40, -1, 1, "Y^{ZNC}_{1}"};
+
+    const AxisSpec axisPsi{18, -PIHalf, PIHalf, "#Psi_{SP}"};
+
+    // Create histograms
+    // Event
+    histos.add("Event/hCent", "FT0C%", kTH1F, {axisCent});
+    histos.add("Event/hVx", "V_{x}", kTH1F, {axisVx});
+    histos.add("Event/hVy", "V_{y}", kTH1F, {axisVy});
+    histos.add("Event/hVz", "V_{z}", kTH1F, {axisVz});
+
+    // Vx / Vy Recent
+    histos.add("Event/hVxVsCent", "<Vx> Vs Cent", kTProfile, {axisCent});
+    histos.add("Event/hVyVsCent", "<Vy> Vs Cent", kTProfile, {axisCent});
+
+    // Gain calib
+    histos.add("QA/GainCalib/hZNASignal", "ZNA Signal", kTH2F, {{4, 0, 4}, {axisZDCEnergy}});
+    histos.add("QA/GainCalib/hZNCSignal", "ZNC Signal", kTH2F, {{4, 0, 4}, {axisZDCEnergy}});
+    histos.add("QA/hZNASignal", "ZNA Signal", kTProfile2D, {{4, 0, 4}, {axisVz}});
+    histos.add("QA/hZNCSignal", "ZNC Signal", kTProfile2D, {{4, 0, 4}, {axisVz}});
+    histos.add("QA/hZNAEnergyCommon", "ZNA Energy Common", kTProfile, {axisVz});
+    histos.add("QA/hZNCEnergyCommon", "ZNC Energy Common", kTProfile, {axisVz});
+
+    // Corrections
+    histos.add("CorrHist/hWtXZNA", "X^{ZNA}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
+    histos.add("CorrHist/hWtYZNA", "Y^{ZNA}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
+    histos.add("CorrHist/hWtXZNC", "X^{ZNC}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
+    histos.add("CorrHist/hWtYZNC", "Y^{ZNC}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
+    histos.add("CorrHist/hUWtXZNA", "X^{ZNA}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
+    histos.add("CorrHist/hUWtYZNA", "Y^{ZNA}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
+    histos.add("CorrHist/hUWtXZNC", "X^{ZNC}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
+    histos.add("CorrHist/hUWtYZNC", "Y^{ZNC}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
+    histos.add("CorrHist/hXZNAVsCent", "X^{ZNA}_{1} Vs Cent", kTProfile, {axisFineCent});
+    histos.add("CorrHist/hXZNAVsVx", "X^{ZNA}_{1} Vs V_{x}", kTProfile, {axisFineVx});
+    histos.add("CorrHist/hXZNAVsVy", "X^{ZNA}_{1} Vs V_{y}", kTProfile, {axisFineVy});
+    histos.add("CorrHist/hXZNAVsVz", "X^{ZNA}_{1} Vs V_{z}", kTProfile, {axisFineVz});
+    histos.add("CorrHist/hYZNAVsCent", "Y^{ZNA}_{1} Vs Cent", kTProfile, {axisFineCent});
+    histos.add("CorrHist/hYZNAVsVx", "Y^{ZNA}_{1} Vs V_{x}", kTProfile, {axisFineVx});
+    histos.add("CorrHist/hYZNAVsVy", "Y^{ZNA}_{1} Vs V_{y}", kTProfile, {axisFineVy});
+    histos.add("CorrHist/hYZNAVsVz", "Y^{ZNA}_{1} Vs V_{z}", kTProfile, {axisFineVz});
+    histos.add("CorrHist/hXZNCVsCent", "X^{ZNC}_{1} Vs Cent", kTProfile, {axisFineCent});
+    histos.add("CorrHist/hXZNCVsVx", "X^{ZNC}_{1} Vs V_{x}", kTProfile, {axisFineVx});
+    histos.add("CorrHist/hXZNCVsVy", "X^{ZNC}_{1} Vs V_{y}", kTProfile, {axisFineVy});
+    histos.add("CorrHist/hXZNCVsVz", "X^{ZNC}_{1} Vs V_{z}", kTProfile, {axisFineVz});
+    histos.add("CorrHist/hYZNCVsCent", "Y^{ZNC}_{1} Vs Cent", kTProfile, {axisFineCent});
+    histos.add("CorrHist/hYZNCVsVx", "Y^{ZNC}_{1} Vs V_{x}", kTProfile, {axisFineVx});
+    histos.add("CorrHist/hYZNCVsVy", "Y^{ZNC}_{1} Vs V_{y}", kTProfile, {axisFineVy});
+    histos.add("CorrHist/hYZNCVsVz", "Y^{ZNC}_{1} Vs V_{z}", kTProfile, {axisFineVz});
+
+    // Checks
+    histos.add("Checks/hPsiSPA", "#Psi_{SP}^{A} distribution", kTH2F, {axisCent, axisPsi});
+    histos.add("Checks/hPsiSPC", "#Psi_{SP}^{C} distribution", kTH2F, {axisCent, axisPsi});
+    histos.add("Checks/hCosPsiSPAC", "Cos(#Psi_{SP}^{A} #minus #Psi_{SP}^{C}) distribution", kTProfile, {axisCent});
+    histos.add("Checks/hSinPsiSPAC", "Sin(#Psi_{SP}^{A} #minus #Psi_{SP}^{C}) distribution", kTProfile, {axisCent});
+    histos.add("Checks/hXaXc", "X^{A}_{1}X^{C}_{1}", kTProfile, {axisCent});
+    histos.add("Checks/hYaYc", "Y^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent});
+    histos.add("Checks/hXaYc", "X^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent});
+    histos.add("Checks/hYaXc", "Y^{A}_{1}X^{C}_{1}", kTProfile, {axisCent});
+
+    // Directed flow QXY vector
+    histos.add("DF/hQaQc", "X^{A}_{1}X^{C}_{1} + Y^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent});
   }
 
+  // Select collsion
   template <typename C>
   bool selCollision(C const& col)
   {
@@ -261,6 +512,293 @@ struct FlowEventPlane {
     return true;
   }
 
+  // Load Gain Calibrations and ZDC Q-Vector Recentering Corrections
+  void loadCorrections()
+  {
+    // Load Vx / Vy recentering
+    if (cRecentVxVy) {
+      std::string ccdbPath = static_cast<std::string>(cCcdbPath) + "/VxVyRecent" + "/Run" + std::to_string(cRunNum);
+      auto ccdbObj = ccdbService->getForTimeStamp<TList>(ccdbPath, nolaterthan.value);
+      CorrectionHistContainer.hVx = reinterpret_cast<TProfile*>(ccdbObj->FindObject("hVx"));
+      CorrectionHistContainer.hVy = reinterpret_cast<TProfile*>(ccdbObj->FindObject("hVy"));
+    }
+
+    // Load ZDC gain calibration
+    if (cDoGainCalib) {
+      std::string ccdbPath = static_cast<std::string>(cCcdbPath) + "/GainCalib" + "/Run" + std::to_string(cRunNum);
+      auto ccdbObj = ccdbService->getForTimeStamp<TList>(ccdbPath, nolaterthan.value);
+      CorrectionHistContainer.hGainCalib[0] = reinterpret_cast<TH2F*>(ccdbObj->FindObject("hZNASignal"));
+      CorrectionHistContainer.hGainCalib[1] = reinterpret_cast<TH2F*>(ccdbObj->FindObject("hZNCSignal"));
+    }
+
+    // Load shift corrections for ZDC Q-Vectors
+    if (cApplyRecentCorr) {
+      std::vector<int> vCorrFlags = static_cast<std::vector<int>>(cCorrFlagVector);
+      int nitr = vCorrFlags.size();
+      CorrectionType corrType = kFineCorr;
+
+      for (int i = 0; i < nitr; ++i) {
+        // Skip correction if corrFlag != 1
+        if (vCorrFlags[i] != 1) {
+          continue;
+        }
+
+        // Set correction type
+        if (i % kNCorr == 0) {
+          corrType = kCoarseCorr;
+        } else {
+          corrType = kFineCorr;
+        }
+
+        // Set ccdb path
+        std::string ccdbPath = static_cast<std::string>(cCcdbPath) + "/CorrItr_" + std::to_string(i + 1) + "/Run" + std::to_string(cRunNum);
+
+        // Get object from CCDB
+        auto ccdbObject = ccdbService->getForTimeStamp<TList>(ccdbPath, nolaterthan.value);
+
+        // Check CCDB Object
+        if (!ccdbObject) {
+          LOGF(warning, "CCDB OBJECT NOT FOUND");
+          return;
+        }
+
+        // Store histograms in Hist Container
+        std::vector<std::vector<std::string>> vHistNames = corrTypeHistNameMap.at(corrType);
+        int cntrx = 0;
+        for (auto const& x : vHistNames) {
+          int cntry = 0;
+          for (auto const& y : x) {
+            if (corrType == kFineCorr) {
+              CorrectionHistContainer.vFineCorrHist[i][cntrx][cntry] = reinterpret_cast<TProfile*>(ccdbObject->FindObject(y.c_str()));
+            } else {
+              CorrectionHistContainer.vCoarseCorrHist[i][cntrx][cntry] = reinterpret_cast<THnSparseF*>(ccdbObject->FindObject(y.c_str()));
+            }
+            ++cntry;
+          }
+          ++cntrx;
+        }
+      }
+    }
+  }
+
+  // Apply gain calibrations
+  void gainCalib(float const& vz, std::array<float, 4>& eA, std::array<float, 4>& eC)
+  {
+    float vA = 0., vC = 0.;
+    for (int i = 0; i < static_cast<int>(eA.size()); ++i) {
+      vA = CorrectionHistContainer.hGainCalib[0]->GetBinContent(CorrectionHistContainer.hGainCalib[0]->FindBin(i + 0.5, vz));
+      vC = CorrectionHistContainer.hGainCalib[1]->GetBinContent(CorrectionHistContainer.hGainCalib[1]->FindBin(i + 0.5, vz));
+      eA[i] *= vA;
+      eC[i] *= vC;
+    }
+  }
+
+  std::vector<float> getAvgCorrFactors(int const& itr, CorrectionType const& corrType, std::array<float, 4> const& vCollParam)
+  {
+    std::vector<float> vAvgOutput = {0., 0., 0., 0.};
+    int binarray[4];
+    if (corrType == kCoarseCorr) {
+      int cntrx = 0;
+      for (auto const& v : CorrectionHistContainer.vCoarseCorrHist[itr]) {
+        for (auto const& h : v) {
+          binarray[kCent] = h->GetAxis(kCent)->FindBin(vCollParam[kCent]);
+          binarray[kVx] = h->GetAxis(kVx)->FindBin(vCollParam[kVx]);
+          binarray[kVy] = h->GetAxis(kVy)->FindBin(vCollParam[kVy]);
+          binarray[kVz] = h->GetAxis(kVz)->FindBin(vCollParam[kVz]);
+          vAvgOutput[cntrx] += h->GetBinContent(h->GetBin(binarray));
+        }
+        ++cntrx;
+      }
+    } else {
+      int cntrx = 0;
+      for (auto const& v : CorrectionHistContainer.vFineCorrHist[itr]) {
+        int cntry = 0;
+        for (auto const& h : v) {
+          vAvgOutput[cntrx] += h->GetBinContent(h->GetXaxis()->FindBin(vCollParam[cntry]));
+          ++cntry;
+        }
+        ++cntrx;
+      }
+    }
+
+    return vAvgOutput;
+  }
+
+  void applyCorrection(const std::array<float, 4>& inputParam, std::array<float, 4>& outputParam)
+  {
+    std::vector<int> vCorrFlags = static_cast<std::vector<int>>(cCorrFlagVector);
+    int nitr = vCorrFlags.size();
+    CorrectionType corrType = kFineCorr;
+
+    // Correction iterations
+    for (int i = 0; i < nitr; ++i) {
+      // Don't correct if corrFlag != 1
+      if (vCorrFlags[i] != 1) {
+        continue;
+      }
+
+      // Set correction type
+      if (i % kNCorr == 0) {
+        corrType = kCoarseCorr;
+      } else {
+        corrType = kFineCorr;
+      }
+
+      // Get averages
+      std::vector<float> vAvg = getAvgCorrFactors(i, corrType, inputParam);
+
+      // Apply recentering
+      outputParam[kXa] -= vAvg[kXa];
+      outputParam[kYa] -= vAvg[kYa];
+      outputParam[kXc] -= vAvg[kXc];
+      outputParam[kYc] -= vAvg[kYc];
+    }
+  }
+
+  void fillCorrHist(std::array<float, 4> const& vCollParam, std::array<float, 4> const& vSP)
+  {
+    histos.fill(HIST("CorrHist/hWtXZNA"), vCollParam[kCent], vCollParam[kVx], vCollParam[kVy], vCollParam[kVz], vSP[kXa]);
+    histos.fill(HIST("CorrHist/hWtYZNA"), vCollParam[kCent], vCollParam[kVx], vCollParam[kVy], vCollParam[kVz], vSP[kYa]);
+    histos.fill(HIST("CorrHist/hWtXZNC"), vCollParam[kCent], vCollParam[kVx], vCollParam[kVy], vCollParam[kVz], vSP[kXc]);
+    histos.fill(HIST("CorrHist/hWtYZNC"), vCollParam[kCent], vCollParam[kVx], vCollParam[kVy], vCollParam[kVz], vSP[kYc]);
+    histos.fill(HIST("CorrHist/hUWtXZNA"), vCollParam[kCent], vCollParam[kVx], vCollParam[kVy], vCollParam[kVz]);
+    histos.fill(HIST("CorrHist/hUWtYZNA"), vCollParam[kCent], vCollParam[kVx], vCollParam[kVy], vCollParam[kVz]);
+    histos.fill(HIST("CorrHist/hUWtXZNC"), vCollParam[kCent], vCollParam[kVx], vCollParam[kVy], vCollParam[kVz]);
+    histos.fill(HIST("CorrHist/hUWtYZNC"), vCollParam[kCent], vCollParam[kVx], vCollParam[kVy], vCollParam[kVz]);
+    histos.fill(HIST("CorrHist/hXZNAVsCent"), vCollParam[kCent], vSP[kXa]);
+    histos.fill(HIST("CorrHist/hXZNAVsVx"), vCollParam[kVx], vSP[kXa]);
+    histos.fill(HIST("CorrHist/hXZNAVsVy"), vCollParam[kVy], vSP[kXa]);
+    histos.fill(HIST("CorrHist/hXZNAVsVz"), vCollParam[kVz], vSP[kXa]);
+    histos.fill(HIST("CorrHist/hYZNAVsCent"), vCollParam[kCent], vSP[kYa]);
+    histos.fill(HIST("CorrHist/hYZNAVsVx"), vCollParam[kVx], vSP[kYa]);
+    histos.fill(HIST("CorrHist/hYZNAVsVy"), vCollParam[kVy], vSP[kYa]);
+    histos.fill(HIST("CorrHist/hYZNAVsVz"), vCollParam[kVz], vSP[kYa]);
+    histos.fill(HIST("CorrHist/hXZNCVsCent"), vCollParam[kCent], vSP[kXc]);
+    histos.fill(HIST("CorrHist/hXZNCVsVx"), vCollParam[kVx], vSP[kXc]);
+    histos.fill(HIST("CorrHist/hXZNCVsVy"), vCollParam[kVy], vSP[kXc]);
+    histos.fill(HIST("CorrHist/hXZNCVsVz"), vCollParam[kVz], vSP[kXc]);
+    histos.fill(HIST("CorrHist/hYZNCVsCent"), vCollParam[kCent], vSP[kYc]);
+    histos.fill(HIST("CorrHist/hYZNCVsVx"), vCollParam[kVx], vSP[kYc]);
+    histos.fill(HIST("CorrHist/hYZNCVsVy"), vCollParam[kVy], vSP[kYc]);
+    histos.fill(HIST("CorrHist/hYZNCVsVz"), vCollParam[kVz], vSP[kYc]);
+  }
+
+  template <typename C, typename B>
+  bool analyzeCollision(B const& bc, C const& collision, std::array<float, 4>& vSP)
+  {
+    // Check ZDC
+    if (!bc.has_zdc()) {
+      return false;
+    }
+
+    // Event selection
+    if (!selCollision(collision)) {
+      return false;
+    }
+    posX = collision.posX();
+    posY = collision.posY();
+    posZ = collision.posZ();
+    std::array<float, 4> vCollParam = {cent, posX, posY, posZ};
+
+    // Zdc information
+    auto zdc = bc.zdc();
+    auto znaEnergy = zdc.energySectorZNA();
+    auto zncEnergy = zdc.energySectorZNC();
+    auto znaEnergyCommon = zdc.energyCommonZNA();
+    auto zncEnergyCommon = zdc.energyCommonZNC();
+
+    // Check energy deposits
+    if (znaEnergyCommon <= 0 || zncEnergyCommon <= 0 || znaEnergy[0] <= 0 || znaEnergy[1] <= 0 || znaEnergy[2] <= 0 || znaEnergy[3] <= 0 || zncEnergy[0] <= 0 || zncEnergy[1] <= 0 || zncEnergy[2] <= 0 || zncEnergy[3] <= 0) {
+      return false;
+    }
+
+    // Selected Collision with required ZDC
+    // Fill avg Vx / Vy
+    histos.fill(HIST("Event/hVxVsCent"), cent, posX);
+    histos.fill(HIST("Event/hVyVsCent"), cent, posY);
+
+    // Apply Vx / Vy recentering
+    if (cRecentVxVy) {
+      vCollParam[kVx] -= CorrectionHistContainer.hVx->GetBinContent(CorrectionHistContainer.hVx->FindBin(cent));
+      vCollParam[kVy] -= CorrectionHistContainer.hVy->GetBinContent(CorrectionHistContainer.hVy->FindBin(cent));
+    }
+
+    // Apply Vx / Vy selection [-10., 10.] mm
+    if (std::abs(vCollParam[kVx]) >= cAxisVxMax || std::abs(vCollParam[kVy]) >= cAxisVyMax) {
+      return false;
+    }
+
+    // Fill gain calib histograms
+    histos.fill(HIST("QA/hZNAEnergyCommon"), vCollParam[kVz], znaEnergyCommon);
+    histos.fill(HIST("QA/hZNCEnergyCommon"), vCollParam[kVz], zncEnergyCommon);
+    for (int iCh = 0; iCh < kXYAC; ++iCh) {
+      histos.fill(HIST("QA/hZNASignal"), iCh + 0.5, vCollParam[kVz], znaEnergy[iCh]);
+      histos.fill(HIST("QA/hZNCSignal"), iCh + 0.5, vCollParam[kVz], zncEnergy[iCh]);
+    }
+
+    // Do gain calibration
+    if (cDoGainCalib) {
+      gainCalib(vCollParam[kVz], znaEnergy, zncEnergy);
+    }
+
+    // Fill zdc signal
+    for (int iCh = 0; iCh < kXYAC; ++iCh) {
+      histos.fill(HIST("QA/GainCalib/hZNASignal"), iCh + 0.5, znaEnergy[iCh]);
+      histos.fill(HIST("QA/GainCalib/hZNCSignal"), iCh + 0.5, zncEnergy[iCh]);
+    }
+
+    // Fill event QA
+    histos.fill(HIST("Event/hCent"), vCollParam[kCent]);
+    histos.fill(HIST("Event/hVx"), vCollParam[kVx]);
+    histos.fill(HIST("Event/hVy"), vCollParam[kVy]);
+    histos.fill(HIST("Event/hVz"), vCollParam[kVz]);
+
+    auto alphaZDC = 0.395;
+    const double x[4] = {-1.75, 1.75, -1.75, 1.75};
+    const double y[4] = {-1.75, -1.75, 1.75, 1.75};
+
+    // Calculate X and Y
+    float znaXNum = 0., znaYNum = 0., zncXNum = 0., zncYNum = 0.;
+    float znaDen = 0., zncDen = 0.;
+    float znaWt = 0., zncWt = 0.;
+
+    // Loop over zdc sectors
+    for (int i = 0; i < kXYAC; ++i) {
+      if (cUseAlphaZDC) {
+        znaWt = std::pow(znaEnergy[i], alphaZDC);
+        zncWt = std::pow(zncEnergy[i], alphaZDC);
+      } else {
+        znaWt = znaEnergy[i];
+        zncWt = zncEnergy[i];
+      }
+      znaXNum -= znaWt * x[i];
+      znaYNum += znaWt * y[i];
+      zncXNum += zncWt * x[i];
+      zncYNum += zncWt * y[i];
+      znaDen += znaWt;
+      zncDen += zncWt;
+    }
+
+    if (znaDen < zdcDenThrs || zncDen < zdcDenThrs) {
+      return false;
+    }
+
+    // Get X and Y for A and C side ZNA
+    vSP[kXa] = znaXNum / znaDen;
+    vSP[kYa] = znaYNum / znaDen;
+    vSP[kXc] = zncXNum / zncDen;
+    vSP[kYc] = zncYNum / zncDen;
+
+    // Do corrections
+    if (cApplyRecentCorr) {
+      applyCorrection(vCollParam, vSP);
+    }
+
+    // Fill X and Y histograms for corrections after each iteration
+    fillCorrHist(vCollParam, vSP);
+    return true;
+  }
+
   // Track Selection
   template <typename T>
   bool selectTrack(T const& track)
@@ -280,77 +818,322 @@ struct FlowEventPlane {
     return true;
   }
 
-  // PID
-  template <typename T>
-  void identifyTrack(T const& track, std::array<float, 4>& nstpc, std::array<float, 4>& nstof)
+  template <ParticleType part1, ParticleType part2, ParticleType part3>
+  bool checkTrackPid(float const& ptCut, float const& trackPt, std::vector<float> const& vTpcNsig, std::vector<float> const& vTofNsig, bool const& tofFlag)
+  {
+    bool retFlag = false;
+    if (tofFlag) {
+      if (vTofNsig[part1] < cTofNSigmaCut && vTofNsig[part2] > cTofRejCut && vTofNsig[part3] > cTofRejCut && vTpcNsig[part1] < cTpcNSigmaCut) {
+        retFlag = true;
+      }
+    } else {
+      if (trackPt < ptCut && vTpcNsig[part1] < cTpcNSigmaCut && vTpcNsig[part2] > cTpcRejCut && vTpcNsig[part3] > cTpcRejCut) {
+        retFlag = true;
+      }
+    }
+    return retFlag;
+  }
+
+  template <ParticleType partType, typename T>
+  bool identifyTrack(T const& track)
   {
     // Electron rejection
     if (std::abs(track.tpcNSigmaPi()) > cTpcRejCut && std::abs(track.tpcNSigmaKa()) > cTpcRejCut && std::abs(track.tpcNSigmaPr()) > cTpcRejCut && track.tpcNSigmaEl() > cTpcElRejCutMin && track.tpcNSigmaEl() < cTpcElRejCutMax) {
-      return;
+      return false;
     }
 
-    // Check track for pion
-    if (std::abs(track.tpcNSigmaPi()) < cTpcNSigmaCut) {
-      nstpc[kPi] = track.tpcNSigmaPi();
-      if (track.hasTOF() && std::abs(track.tofNSigmaPi()) < cTofNSigmaCut) {
-        nstof[kPi] = track.tofNSigmaPi();
-      }
+    // Pion, Kaon, Proton Identification
+    std::vector<float> vPtCut = {cPionPtCut, cKaonPtCut, cProtonPtCut};
+    std::vector<float> vTpcNsig = {std::abs(track.tpcNSigmaPi()), std::abs(track.tpcNSigmaKa()), std::abs(track.tpcNSigmaPr())};
+    std::vector<float> vTofNsig = {std::abs(track.tofNSigmaPi()), std::abs(track.tofNSigmaKa()), std::abs(track.tofNSigmaPr())};
+    bool retFlag = false;
+
+    if (partType == kPi && checkTrackPid<kPi, kKa, kPr>(vPtCut[kPi], track.pt(), vTpcNsig, vTofNsig, track.hasTOF())) {
+      retFlag = true;
+    } else if (partType == kKa && checkTrackPid<kKa, kPi, kPr>(vPtCut[kKa], track.pt(), vTpcNsig, vTofNsig, track.hasTOF())) {
+      retFlag = true;
+    } else if (partType == kPr && checkTrackPid<kPr, kPi, kKa>(vPtCut[kPr], track.pt(), vTpcNsig, vTofNsig, track.hasTOF())) {
+      retFlag = true;
+    } else {
+      return false;
     }
 
-    // Check track for kaon
-    if (std::abs(track.tpcNSigmaKa()) < cTpcNSigmaCut) {
-      nstpc[kKa] = track.tpcNSigmaKa();
-      if (track.hasTOF() && std::abs(track.tofNSigmaKa()) < cTofNSigmaCut) {
-        nstof[kKa] = track.tofNSigmaKa();
-      }
+    return retFlag;
+  }
+
+  using BCsRun3 = soa::Join<aod::BCsWithTimestamps, aod::Run3MatchedToBCSparse>;
+  using CollisionsRun3 = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::CentFT0Ms, aod::CentFV0As, aod::TPCMults, aod::PVMults, aod::MultsGlobal, aod::MultsExtra>;
+  using Tracks = soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::TOFSignal, aod::pidTPCEl, aod::pidTPCPi, aod::pidTOFPi, aod::pidTPCKa, aod::pidTOFKa, aod::pidTPCPr, aod::pidTOFPr>;
+
+  void processDummy(CollisionsRun3::iterator const&) {}
+
+  PROCESS_SWITCH(SpectatorPlaneTableProducer, processDummy, "Dummy process", true);
+
+  void processSpectatorPlane(CollisionsRun3::iterator const& collision, BCsRun3 const&, aod::Zdcs const&)
+  {
+    // Get bunch crossing
+    auto bc = collision.template foundBC_as<BCsRun3>();
+    cRunNum = collision.template foundBC_as<BCsRun3>().runNumber();
+
+    if (cLoadCorrection && lRunNum != cRunNum) {
+      loadCorrections();
     }
 
-    // Check track for proton
-    if (std::abs(track.tpcNSigmaPr()) < cTpcNSigmaCut) {
-      nstpc[kPr] = track.tpcNSigmaPr();
-      if (track.hasTOF() && std::abs(track.tofNSigmaPr()) < cTofNSigmaCut) {
-        nstof[kPr] = track.tofNSigmaPr();
-      }
+    // Analyze collision and get Spectator Plane Vector
+    std::array<float, 4> vSP = {0., 0., 0., 0.};
+    bool colSPExtFlag = analyzeCollision(bc, collision, vSP);
+
+    // Update run number
+    lRunNum = cRunNum;
+
+    // Fill histograms if SP flag is true
+    if (colSPExtFlag) {
+      // Evaluate spectator plane angle and [X,Y] correlations
+      float psiA = std::atan2(vSP[kYa], vSP[kXa]);
+      float psiC = std::atan2(vSP[kYc], vSP[kXc]);
+      histos.fill(HIST("Checks/hPsiSPA"), cent, psiA);
+      histos.fill(HIST("Checks/hPsiSPC"), cent, psiC);
+      histos.fill(HIST("Checks/hCosPsiSPAC"), cent, std::cos(psiA - psiC));
+      histos.fill(HIST("Checks/hSinPsiSPAC"), cent, std::sin(psiA - psiC));
+      histos.fill(HIST("Checks/hXaXc"), cent, (vSP[kXa] * vSP[kXc]));
+      histos.fill(HIST("Checks/hYaYc"), cent, (vSP[kYa] * vSP[kYc]));
+      histos.fill(HIST("Checks/hXaYc"), cent, (vSP[kXa] * vSP[kYc]));
+      histos.fill(HIST("Checks/hYaXc"), cent, (vSP[kYa] * vSP[kXc]));
+
+      // Directed flow QXY vector
+      float qac = (vSP[kXa] * vSP[kXc]) + (vSP[kYa] * vSP[kYc]);
+      histos.fill(HIST("DF/hQaQc"), cent, qac);
     }
+
+    // Fill table
+    colSPExtTable(colSPExtFlag, vSP[kXa], vSP[kYa], vSP[kXc], vSP[kYc]);
   }
 
-  // Id track table
-  template <typename T>
-  void analyzeIdHadrons(T const& tracks)
+  PROCESS_SWITCH(SpectatorPlaneTableProducer, processSpectatorPlane, "Spectator Plane Process", false);
+
+  void processIdHadrons(Tracks const& tracks)
   {
     for (auto const& track : tracks) {
-      // Global track selection
-      if (!selectTrack(track)) {
-        continue;
-      }
+      bool chargedFlag = selectTrack(track);
+      bool pionFlag = identifyTrack<kPi>(track);
+      bool kaonFlag = identifyTrack<kKa>(track);
+      bool protonFlag = identifyTrack<kPr>(track);
+      tracksIdTable(chargedFlag, pionFlag, kaonFlag, protonFlag);
+    }
+  }
+
+  PROCESS_SWITCH(SpectatorPlaneTableProducer, processIdHadrons, "Hadron Id Process", false);
+};
+
+struct FlowEventPlane {
+  // Tracks
+  Configurable<int> cEtaBins{"cEtaBins", 5, "# of eta bins"};
+  Configurable<float> cEtaCut{"cEtaCut", 0.8, "Rapidity cut"};
+
+  // Pi,Ka,Pr
+  Configurable<float> cMinPtPi{"cMinPtPi", 0.2, "Pion min pT"};
+  Configurable<float> cMinPtKa{"cMinPtKa", 0.3, "Kaon min pT"};
+  Configurable<float> cMinPtPr{"cMinPtPr", 0.5, "Proton min pT"};
+
+  // Resonance
+  Configurable<int> cPhiInvMassBins{"cPhiInvMassBins", 500, "# of Phi mass bins"};
+
+  // V0
+  // Tracks
+  Configurable<float> cTrackPtCut{"cTrackPtCut", 0.1, "p_{T} minimum"};
+  Configurable<float> cTrackEtaCut{"cTrackEtaCut", 0.8, "Pseudorapidity cut"};
+  Configurable<double> cTpcNsigmaCut{"cTpcNsigmaCut", 3.0, "TPC NSigma Selection Cut"};
+
+  // V0s
+  Configurable<int> cV0TypeSelection{"cV0TypeSelection", 1, "V0 Type Selection"};
+  Configurable<float> cMinDcaProtonToPV{"cMinDcaProtonToPV", 0.01, "Minimum Proton DCAr to PV"};
+  Configurable<float> cMinDcaPionToPV{"cMinDcaPionToPV", 0.1, "Minimum Pion DCAr to PV"};
+  Configurable<float> cDcaV0Dau{"cDcaV0Dau", 1., "DCA between V0 daughters"};
+  Configurable<float> cMinV0Radius{"cMinV0Radius", 0.5, "Minimum V0 radius from PV"};
+  Configurable<float> cK0ShortCTau{"cK0ShortCTau", 20.0, "Decay length cut K0Short"};
+  Configurable<float> cLambdaCTau{"cLambdaCTau", 30.0, "Decay length cut Lambda"};
+  Configurable<float> cK0ShortCosPA{"cK0ShortCosPA", 0.97, "K0Short CosPA"};
+  Configurable<float> cLambdaCosPA{"cLambdaCosPA", 0.98, "Lambda CosPA"};
+  Configurable<float> cK0SMassRej{"cK0SMassRej", 0.01, "Reject K0Short Candidates"};
+  Configurable<float> cArmPodSel{"cArmPodSel", 0.2, "Armentros-Podolanski Selection for K0S"};
+  Configurable<float> cK0SMinPt{"cK0SMinPt", 0.4, "K0S Min pT"};
+  Configurable<float> cK0SMaxPt{"cK0SMaxPt", 6.0, "K0S Max pT"};
+  Configurable<float> cLambdaMinPt{"cLambdaMinPt", 0.6, "Lambda Min pT"};
+  Configurable<float> cLambdaMaxPt{"cLambdaMaxPt", 6.0, "Lambda Max pT"};
 
-      // Dca
-      std::array<float, 2> dca = {track.dcaXY(), track.dcaZ()};
+  // Histogram registry: an object to hold your histograms
+  HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
+
+  // Global objects
+  float cent = 0.;
+  std::array<float, 4> vSP = {0., 0., 0., 0.};
+  std::map<ResoType, std::array<float, 2>> mResoDauMass = {{kPhi0, {MassKaonCharged, MassKaonCharged}}, {kKStar, {MassPionCharged, MassKaonCharged}}};
+  std::map<ResoType, float> mResoMass = {{kPhi0, MassPhi}, {kKStar, MassKaonCharged}};
+
+  void init(InitContext const&)
+  {
+    // Define axes
+    const AxisSpec axisCent{100, 0., 100, "FT0C%"};
+
+    const AxisSpec axisXYac{600, -6, 6, "Q^{t}Q^{p}"};
+    const AxisSpec axisV1{400, -4, 4, "v_{1}"};
+
+    const AxisSpec axisTrackPt{100, 0., 10., "p_{T} (GeV/#it{c})"};
+    const AxisSpec axisTrackEta{cEtaBins, -0.8, 0.8, "#eta"};
+    const AxisSpec axisTrackDcaXY{60, -0.15, 0.15, "DCA_{XY}"};
+    const AxisSpec axisTrackDcaZ{230, -1.15, 1.15, "DCA_{XY}"};
+    const AxisSpec axisTrackdEdx{360, 20, 200, "#frac{dE}{dx}"};
+    const AxisSpec axisTrackNSigma{161, -4.025, 4.025, {"n#sigma"}};
+
+    const AxisSpec axisPhiInvMass{cPhiInvMassBins, 0.99, 1.12, "M_{KK} (GeV/#it{c}^{2}"};
+    const AxisSpec axisMomPID(80, 0, 4, "p_{T} (GeV/#it{c})");
+    const AxisSpec axisNsigma(401, -10.025, 10.025, {"n#sigma"});
+    const AxisSpec axisdEdx(360, 20, 200, "#frac{dE}{dx}");
+    const AxisSpec axisTrackTofSignal(240, 0, 1.2, "#beta");
+    const AxisSpec axisRadius(2000, 0, 200, "r(cm)");
+    const AxisSpec axisCosPA(300, 0.97, 1.0, "cos(#theta_{PA})");
+    const AxisSpec axisDcaV0PV(1000, 0., 10., "dca (cm)");
+    const AxisSpec axisDcaProngPV(5000, -50., 50., "dca (cm)");
+    const AxisSpec axisDcaDau(75, 0., 1.5, "Daug DCA (#sigma)");
+    const AxisSpec axisCTau(2000, 0, 200, "c#tau (cm)");
+    const AxisSpec axisAlpha(40, -1, 1, "#alpha");
+    const AxisSpec axisQtarm(40, 0, 0.4, "q_{T}");
+    const AxisSpec axisLambdaInvMass(140, 1.08, 1.15, "M_{p#pi} (GeV/#it{c}^{2})");
+    const AxisSpec axisK0ShortInvMass(200, 0.4, 0.6, "M_{#pi#pi} (GeV/#it{c}^{2})");
+
+    // Create histograms
+    // Charged particles
+    if (doprocessChargedFlow) {
+      histos.add("TrackQA/hPtDcaXY", "DCA_{XY} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaXY});
+      histos.add("TrackQA/hPtDcaZ", "DCA_{Z} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaZ});
+      histos.add("TrackQA/hTrackTPCdEdX", "hTrackTPCdEdX", kTH2F, {axisMomPID, axisdEdx});
+      histos.add("DF/hAQu", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta});
+      histos.add("DF/hCQu", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta});
+      histos.add("DF/hAQuPos", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta});
+      histos.add("DF/hCQuPos", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta});
+      histos.add("DF/hAQuNeg", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta});
+      histos.add("DF/hCQuNeg", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta});
+    }
+
+    // Identified hadrons
+    if (doprocessIdFlow) {
+      histos.add("PartId/Pion/hdEdX", "dE/dx vs pT", kTH2F, {axisMomPID, axisTrackdEdx});
+      histos.add("PartId/Pion/hTOFSignal", "#beta_{TOF} vs p_{T}", kTH2F, {axisMomPID, axisTrackTofSignal});
+      histos.add("PartId/Pion/hTPCNSigma", "n#sigma_{TPC} vs p_{T}", kTH2F, {axisMomPID, axisTrackNSigma});
+      histos.add("PartId/Pion/hTOFNSigma", "n#sigma_{TOF} vs p_{T}", kTH2F, {axisMomPID, axisTrackNSigma});
+      histos.add("PartId/Pion/hAQuPos", "PartId/Pion/hAQuPos", kTProfile2D, {axisCent, axisTrackEta});
+      histos.add("PartId/Pion/hAQuNeg", "PartId/Pion/hAQuNeg", kTProfile2D, {axisCent, axisTrackEta});
+      histos.add("PartId/Pion/hCQuPos", "PartId/Pion/hCQuPos", kTProfile2D, {axisCent, axisTrackEta});
+      histos.add("PartId/Pion/hCQuNeg", "PartId/Pion/hCQuNeg", kTProfile2D, {axisCent, axisTrackEta});
+      histos.addClone("PartId/Pion/", "PartId/Kaon/");
+      histos.addClone("PartId/Pion/", "PartId/Proton/");
+    }
+
+    // Resonance
+    if (doprocessResoFlow) {
+      histos.add("Reso/Phi/hSigCentEtaInvMass", "hUSCentEtaInvMass", kTH3F, {axisCent, axisTrackEta, axisPhiInvMass});
+      histos.add("Reso/Phi/hBkgCentEtaInvMass", "hLSCentEtaInvMass", kTH3F, {axisCent, axisTrackEta, axisPhiInvMass});
+      histos.add("Reso/Phi/Sig/hQuA", "hPhiQuA", kTProfile3D, {axisCent, axisTrackEta, axisPhiInvMass});
+      histos.add("Reso/Phi/Sig/hQuC", "hPhiQuC", kTProfile3D, {axisCent, axisTrackEta, axisPhiInvMass});
+      histos.add("Reso/Phi/Bkg/hQuA", "hPhiQuA", kTProfile3D, {axisCent, axisTrackEta, axisPhiInvMass});
+      histos.add("Reso/Phi/Bkg/hQuC", "hPhiQuC", kTProfile3D, {axisCent, axisTrackEta, axisPhiInvMass});
+    }
+
+    // Lambda
+    if (doprocessV0Flow) {
+      histos.add("V0/Lambda/QA/hQtVsAlpha", "Armentros-Podolanski Plot", kTH2F, {axisAlpha, axisQtarm});
+      histos.add("V0/Lambda/QA/hDcaV0Dau", "DCA between V0 daughters", kTH1F, {axisDcaDau});
+      histos.add("V0/Lambda/QA/hDcaPosToPv", "DCA positive prong to PV", kTH1F, {axisDcaProngPV});
+      histos.add("V0/Lambda/QA/hDcaNegToPv", "DCA negative prong to PV", kTH1F, {axisDcaProngPV});
+      histos.add("V0/Lambda/QA/hDcaV0ToPv", "DCA V0 to PV", kTH1F, {axisDcaV0PV});
+      histos.add("V0/Lambda/QA/hCosPa", "cos(#theta_{PA})", kTH1F, {axisCosPA});
+      histos.add("V0/Lambda/QA/hRxy", "V_{0} Decay Radius in XY plane", kTH1F, {axisRadius});
+      histos.add("V0/Lambda/QA/hCTau", "V_{0} c#tau", kTH1F, {axisCTau});
+      histos.add("V0/Lambda/QA/hPosdEdXVsP", "TPC Signal Pos-Prong", kTH2F, {axisMomPID, axisdEdx});
+      histos.add("V0/Lambda/QA/hNegdEdXVsP", "TPC Signal Neg-Prong", kTH2F, {axisMomPID, axisdEdx});
+      histos.add("V0/Lambda/QA/hPosNsigPrVsP", "TPC n#sigma Pos Prong", kTH2F, {axisMomPID, axisNsigma});
+      histos.add("V0/Lambda/QA/hNegNsigPrVsP", "TPC n#sigma Neg Prong", kTH2F, {axisMomPID, axisNsigma});
+      histos.add("V0/Lambda/QA/hPosNsigPiVsP", "TPC n#sigma Pos Prong", kTH2F, {axisMomPID, axisNsigma});
+      histos.add("V0/Lambda/QA/hNegNsigPiVsP", "TPC n#sigma Neg Prong", kTH2F, {axisMomPID, axisNsigma});
+      histos.addClone("V0/Lambda/", "V0/K0Short/");
+      histos.add("V0/Lambda/hMassVsRap", "hMassVsRap", kTH3F, {axisCent, axisTrackEta, axisLambdaInvMass});
+      histos.add("V0/Lambda/Flow/hQuA", "hQuA", kTProfile3D, {axisCent, axisTrackEta, axisLambdaInvMass});
+      histos.add("V0/Lambda/Flow/hQuC", "hQuC", kTProfile3D, {axisCent, axisTrackEta, axisLambdaInvMass});
+      histos.addClone("V0/Lambda/", "V0/AntiLambda/");
+      histos.addClone("V0/Lambda/", "V0/LambdaAntiLambda/");
+      histos.add("V0/K0Short/hMassVsRap", "hMassVsRap", kTH3F, {axisCent, axisTrackEta, axisK0ShortInvMass});
+      histos.add("V0/K0Short/Flow/hQuA", "hQuA", kTProfile3D, {axisCent, axisTrackEta, axisK0ShortInvMass});
+      histos.add("V0/K0Short/Flow/hQuC", "hQuC", kTProfile3D, {axisCent, axisTrackEta, axisK0ShortInvMass});
+    }
+  }
+
+  template <typename C>
+  bool selCollision(C const& collision, std::array<float, 4>& v)
+  {
+    // Check collision
+    if (!collision.selColFlag()) {
+      return false;
+    }
 
-      // QA plots
-      histos.fill(HIST("QA/hDcaXY"), track.pt(), track.dcaXY());
-      histos.fill(HIST("QA/hDcaZ"), track.pt(), track.dcaZ());
-      histos.fill(HIST("QA/hdEdX"), track.pt(), track.tpcSignal());
+    // Set centrality
+    cent = collision.centFT0C();
+
+    // Flow vectors
+    v[kXa] = collision.xa();
+    v[kYa] = collision.ya();
+    v[kXc] = collision.xc();
+    v[kYc] = collision.yc();
+
+    return true;
+  }
+
+  template <ParticleType part, typename T>
+  void getIdFlow(T const& tracks)
+  {
+    float ux = 0., uy = 0., v1a = 0., v1c = 0.;
+    float tpcNsigma = 0., tofNsigma = 0.;
+    for (auto const& track : tracks) {
+      static constexpr std::string_view SubDir[] = {"Pion/", "Kaon/", "Proton/"};
+      if (part == kPi && track.pt() > cMinPtPi) {
+        tpcNsigma = track.tpcNSigmaPi();
+        tofNsigma = track.tofNSigmaPi();
+      } else if (part == kKa && track.pt() > cMinPtKa) {
+        tpcNsigma = track.tpcNSigmaKa();
+        tofNsigma = track.tofNSigmaKa();
+      } else if (part == kPr && track.pt() > cMinPtPr) {
+        tpcNsigma = track.tpcNSigmaPr();
+        tofNsigma = track.tofNSigmaPr();
+      } else {
+        return;
+      }
+      histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hdEdX"), track.pt(), track.tpcSignal());
+      histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hTPCNSigma"), track.pt(), tpcNsigma);
       if (track.hasTOF()) {
-        histos.fill(HIST("QA/hTOFSignal"), track.pt(), track.beta());
+        histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hTOFSignal"), track.pt(), track.beta());
+        histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hTOFNSigma"), track.pt(), tofNsigma);
       }
 
-      // Identify track
-      std::array<float, 4> tpc = {track.tpcSignal(), -99., -99., -99.};
-      std::array<float, 4> tof = {track.beta(), -99., -99., -99.};
-      identifyTrack(track, tpc, tof);
+      // Directed flow
+      ux = std::cos(track.phi());
+      uy = std::sin(track.phi());
+      v1a = ux * vSP[kXa] + uy * vSP[kYa];
+      v1c = ux * vSP[kXc] + uy * vSP[kYc];
 
-      // Fill track table
-      tracksIdTable(colIdx, track.px(), track.py(), track.pz(), track.sign(), tpc.data(), tof.data(), dca.data());
+      if (track.sign() > 0) {
+        histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hAQuPos"), cent, track.eta(), v1a);
+        histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hCQuPos"), cent, track.eta(), v1c);
+      } else {
+        histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hAQuNeg"), cent, track.eta(), v1a);
+        histos.fill(HIST("PartId/") + HIST(SubDir[part]) + HIST("hCQuNeg"), cent, track.eta(), v1c);
+      }
     }
   }
 
-  // V0 daughter track selection
   template <V0Type part, typename V, typename T>
-  bool selV0DauTracks(V const& v0, T const& postrack, T const& negtrack, float& tpcNSigmaPosDau, float& tpcNSigmaNegDau)
+  bool selV0DauTracks(V const& v0, T const& postrack, T const& negtrack)
   {
     // Kinematic selection
-    if (postrack.pt() <= cTrackMinPt || negtrack.pt() <= cTrackMinPt || std::abs(postrack.eta()) >= cTrackEtaCut || std::abs(negtrack.eta()) >= cTrackEtaCut) {
+    if (postrack.pt() <= cTrackPtCut || negtrack.pt() <= cTrackPtCut || std::abs(postrack.eta()) >= cTrackEtaCut || std::abs(negtrack.eta()) >= cTrackEtaCut) {
       return false;
     }
 
@@ -371,181 +1154,273 @@ struct FlowEventPlane {
       return false;
     }
 
-    // Daughter track PID
-    switch (part) {
-      // PosDau = Pion, NegDau = Pion
-      case kK0S:
-        tpcNSigmaPosDau = postrack.tpcNSigmaPi();
-        tpcNSigmaNegDau = negtrack.tpcNSigmaPi();
-        break;
+    // Daughter track PID [Dau1 = PosTrack, Dau2 = NegTrack]
+    float tpcNSigmaDau1 = 0., tpcNSigmaDau2 = 0.;
 
-      // PosDau = Proton, NegDau = Pion
+    switch (part) {
+      // Dau1 = Proton, Dau2 = Pion
       case kLambda:
-        tpcNSigmaPosDau = postrack.tpcNSigmaPr();
-        tpcNSigmaNegDau = negtrack.tpcNSigmaPi();
+        tpcNSigmaDau1 = postrack.tpcNSigmaPr();
+        tpcNSigmaDau2 = negtrack.tpcNSigmaPi();
         break;
 
-      // PosDau = Pion, NegDau = Proton
+      // Dau1 = Pion, Dau2 = Proton
       case kAntiLambda:
-        tpcNSigmaPosDau = postrack.tpcNSigmaPi();
-        tpcNSigmaNegDau = negtrack.tpcNSigmaPr();
+        tpcNSigmaDau1 = postrack.tpcNSigmaPi();
+        tpcNSigmaDau2 = negtrack.tpcNSigmaPr();
+        break;
+
+      // Dau1 = Pion, Dau2 = Pion
+      case kK0S:
+        tpcNSigmaDau1 = postrack.tpcNSigmaPi();
+        tpcNSigmaDau2 = negtrack.tpcNSigmaPi();
         break;
     }
 
-    if (std::abs(tpcNSigmaPosDau) >= cTpcNSigmaCut || std::abs(tpcNSigmaNegDau) >= cTpcNSigmaCut) {
+    if (std::abs(tpcNSigmaDau1) >= cTpcNsigmaCut || std::abs(tpcNSigmaDau2) >= cTpcNsigmaCut) {
       return false;
     }
 
-    // All checks passed
     return true;
   }
 
-  // V0 table
-  template <V0Type v0type, typename V, typename T>
-  void fillV0Table(V const& v0, T const&)
+  template <V0Type part, typename C, typename V, typename T>
+  void fillV0QAHist(C const& col, V const& v0, T const&)
   {
-    // V0 parameters
-    float mass = 0., ctau = 0., tpcNSigmaPosDau = 0., tpcNSigmaNegDau = 0.;
+    static constexpr std::string_view SubDir[] = {"V0/K0Short/QA/", "V0/Lambda/QA/", "V0/AntiLambda/QA/"};
+
+    // daugthers
     auto postrack = v0.template posTrack_as<T>();
     auto negtrack = v0.template negTrack_as<T>();
 
-    if (v0type == kK0S) {
-      if (!selV0DauTracks<kK0S>(v0, postrack, negtrack, tpcNSigmaPosDau, tpcNSigmaNegDau)) {
-        return;
-      }
-      mass = v0.mK0Short();
-      ctau = v0.distovertotmom(posX, posY, posZ) * MassKaonNeutral;
-    } else if (v0type == kLambda) {
-      if (!selV0DauTracks<kLambda>(v0, postrack, negtrack, tpcNSigmaPosDau, tpcNSigmaNegDau)) {
-        return;
-      }
-      mass = v0.mLambda();
-      ctau = v0.distovertotmom(posX, posY, posZ) * MassLambda0;
-    } else if (v0type == kAntiLambda) {
-      if (!selV0DauTracks<kAntiLambda>(v0, postrack, negtrack, tpcNSigmaPosDau, tpcNSigmaNegDau)) {
-        return;
-      }
-      mass = v0.mAntiLambda();
-      ctau = v0.distovertotmom(posX, posY, posZ) * MassLambda0;
+    // ctau
+    float mPDG = 0, ctau = 0;
+    if (part == kK0S) {
+      mPDG = MassKaonNeutral;
     } else {
-      return;
+      mPDG = MassLambda0;
     }
-
-    // Apply V0 selection
-    if (ctau > cV0CTau || v0.v0cosPA() < cV0CosPA) {
-      return;
-    }
-
-    // Dca daughters
-    std::array<float, 2> dcaDau = {v0.dcapostopv(), v0.dcanegtopv()};
-
-    // Tpc daughters
-    std::array<float, 2> tpcDau = {tpcNSigmaPosDau, tpcNSigmaNegDau};
-
-    // Fill QA Plots
-    if (v0type == kK0S) {
-      histos.fill(HIST("QA/h2f_armpod_K0s"), v0.alpha(), v0.qtarm());
-    } else if (v0type == kLambda) {
-      histos.fill(HIST("QA/h2f_armpod_Lambda"), v0.alpha(), v0.qtarm());
-    } else if (v0type == kAntiLambda) {
-      histos.fill(HIST("QA/h2f_armpod_AntiLambda"), v0.alpha(), v0.qtarm());
-    }
-
-    // Fill V0 table
-    v0sTable(colIdx, v0.px(), v0.py(), v0.pz(), v0.v0cosPA(), ctau, mass, (uint8_t)v0type, dcaDau.data(), tpcDau.data());
+    ctau = v0.distovertotmom(col.posX(), col.posY(), col.posZ()) * mPDG;
+
+    histos.fill(HIST(SubDir[part]) + HIST("hQtVsAlpha"), v0.alpha(), v0.qtarm());
+    histos.fill(HIST(SubDir[part]) + HIST("hDcaV0Dau"), v0.dcaV0daughters());
+    histos.fill(HIST(SubDir[part]) + HIST("hDcaPosToPv"), v0.dcapostopv());
+    histos.fill(HIST(SubDir[part]) + HIST("hDcaNegToPv"), v0.dcanegtopv());
+    histos.fill(HIST(SubDir[part]) + HIST("hDcaV0ToPv"), v0.dcav0topv());
+    histos.fill(HIST(SubDir[part]) + HIST("hCosPa"), v0.v0cosPA());
+    histos.fill(HIST(SubDir[part]) + HIST("hRxy"), v0.v0radius());
+    histos.fill(HIST(SubDir[part]) + HIST("hCTau"), ctau);
+    histos.fill(HIST(SubDir[part]) + HIST("hPosdEdXVsP"), postrack.tpcInnerParam(), postrack.tpcSignal());
+    histos.fill(HIST(SubDir[part]) + HIST("hNegdEdXVsP"), negtrack.tpcInnerParam(), negtrack.tpcSignal());
+    histos.fill(HIST(SubDir[part]) + HIST("hPosNsigPrVsP"), postrack.tpcInnerParam(), postrack.tpcNSigmaPr());
+    histos.fill(HIST(SubDir[part]) + HIST("hNegNsigPrVsP"), negtrack.tpcInnerParam(), negtrack.tpcNSigmaPr());
+    histos.fill(HIST(SubDir[part]) + HIST("hPosNsigPiVsP"), postrack.tpcInnerParam(), postrack.tpcNSigmaPi());
+    histos.fill(HIST(SubDir[part]) + HIST("hNegNsigPiVsP"), negtrack.tpcInnerParam(), negtrack.tpcNSigmaPi());
   }
 
-  // V0s
-  template <typename V, typename T>
-  void analyzeV0s(V const& V0s, T const& tracks)
+  template <ResoType rt, typename T>
+  void getResoFlow(T const& tracks1, T const& tracks2, std::array<float, 4> const& vSP)
   {
-    for (auto const& v0 : V0s) {
-      // Topological and kinematic selections
-      if (std::abs(v0.eta()) >= cV0EtaCut || v0.dcaV0daughters() >= cDcaV0Dau || v0.v0radius() <= cMinV0Radius || v0.v0Type() != cV0TypeSelection) {
+    float ux = 0., uy = 0., v1a = 0., v1c = 0.;
+    std::array<float, 2> vMassDau = mResoDauMass.at(rt);
+    for (auto const& [track1, track2] : soa::combinations(soa::CombinationsFullIndexPolicy(tracks1, tracks2))) {
+      // Discard same track
+      if (track1.index() == track2.index()) {
         continue;
       }
 
-      // Armenteros-Podolanski Selections
-      histos.fill(HIST("QA/h2f_armpod_before_sel"), v0.alpha(), v0.qtarm());
+      // Discard same charge track
+      if (track1.sign() == track2.sign()) {
+        continue;
+      }
 
-      // K0s, Lambda, Anti-Lambda
-      if (v0.qtarm() >= cArmPodSel * std::abs(v0.alpha())) { // K0s
-        fillV0Table<kK0S>(v0, tracks);
-      } else if ((v0.qtarm() < cArmPodSel * std::abs(v0.alpha())) && (v0.alpha() > 0)) { // Lambda
-        fillV0Table<kLambda>(v0, tracks);
-      } else if ((v0.qtarm() < cArmPodSel * std::abs(v0.alpha())) && (v0.alpha() < 0)) { // Anti-Lambda
-        fillV0Table<kAntiLambda>(v0, tracks);
-      } else {
-        return;
+      // Apply pseudo-rapidity acceptance
+      std::array<float, 3> v = {track1.px() + track2.px(), track1.py() + track2.py(), track1.pz() + track2.pz()};
+      if (std::abs(RecoDecay::eta(v)) >= cEtaCut) {
+        continue;
       }
+
+      // Reconstruct resonance
+      float p = RecoDecay::p((track1.px() + track2.px()), (track1.py() + track2.py()), (track1.pz() + track2.pz()));
+      float e = RecoDecay::e(track1.px(), track1.py(), track1.pz(), vMassDau[0]) + RecoDecay::e(track2.px(), track2.py(), track2.pz(), vMassDau[1]);
+      float m = std::sqrt(RecoDecay::m2(p, e));
+
+      // Get directed flow
+      ux = std::cos(RecoDecay::phi(v));
+      uy = std::sin(RecoDecay::phi(v));
+      v1a = ux * vSP[kXa] + uy * vSP[kYa];
+      v1c = ux * vSP[kXc] + uy * vSP[kYc];
+
+      // Histograms
+      static constexpr std::string_view SubDir[] = {"Reso/Phi/", "Reso/KStar/"};
+
+      // Fill signal histogram
+      histos.fill(HIST(SubDir[rt]) + HIST("hSigCentEtaInvMass"), cent, RecoDecay::eta(v), m);
+      histos.fill(HIST(SubDir[rt]) + HIST("Sig/hQuA"), cent, RecoDecay::eta(v), m, v1a);
+      histos.fill(HIST(SubDir[rt]) + HIST("Sig/hQuC"), cent, RecoDecay::eta(v), m, v1c);
+
+      // Get background
+      p = RecoDecay::p((track1.px() - track2.px()), (track1.py() - track2.py()), (track1.pz() - track2.pz()));
+      m = std::sqrt(RecoDecay::m2(p, e));
+      v[0] = track1.px() - track2.px();
+      v[1] = track1.py() - track2.py();
+      v[2] = track1.pz() - track2.pz();
+      ux = std::cos(RecoDecay::phi(v));
+      uy = std::sin(RecoDecay::phi(v));
+      v1a = ux * vSP[kXa] + uy * vSP[kYa];
+      v1c = ux * vSP[kXc] + uy * vSP[kYc];
+
+      // Fill bkg histogram
+      histos.fill(HIST(SubDir[rt]) + HIST("hBkgCentEtaInvMass"), cent, RecoDecay::eta(v), m);
+      histos.fill(HIST(SubDir[rt]) + HIST("Bkg/hQuA"), cent, RecoDecay::eta(v), m, v1a);
+      histos.fill(HIST(SubDir[rt]) + HIST("Bkg/hQuC"), cent, RecoDecay::eta(v), m, v1c);
     }
   }
 
-  template <typename B, typename C, typename T, typename V>
-  void analyzeCollision(B const& bc, C const& collision, T const& tracks, V const& V0s)
+  using CollisionsRun3 = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::CentFT0Ms, aod::CentFV0As, aod::MultsExtra, aod::ColSPExt>;
+  using Tracks = soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::TOFSignal, aod::pidTOFbeta, aod::pidTPCEl, aod::pidTPCPi, aod::pidTOFPi, aod::pidTPCKa, aod::pidTOFKa, aod::pidTPCPr, aod::pidTOFPr, aod::TracksId>;
+  using TracksV0s = soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::pidTPCPi, aod::pidTPCPr>;
+
+  // Partitions
+  SliceCache cache;
+  Partition<Tracks> chargedTrackPartition = (aod::tracksid::isCharged == true);
+  Partition<Tracks> pionTrackPartition = (aod::tracksid::isCharged == true) && (aod::tracksid::isPion == true);
+  Partition<Tracks> kaonTrackPartition = (aod::tracksid::isCharged == true) && (aod::tracksid::isKaon == true);
+  Partition<Tracks> protonTrackPartition = (aod::tracksid::isCharged == true) && (aod::tracksid::isProton == true);
+
+  void processDummy(aod::Collisions const&) {}
+
+  PROCESS_SWITCH(FlowEventPlane, processDummy, "Dummy process", true);
+
+  void processChargedFlow(CollisionsRun3::iterator const& collision, Tracks const&)
   {
-    // Event selection
-    if (!selCollision(collision)) {
+    // Check collision
+    if (!selCollision(collision, vSP)) {
       return;
     }
-    posX = collision.posX();
-    posY = collision.posY();
-    posZ = collision.posZ();
-
-    // check zdc
-    if (!bc.has_zdc()) {
-      return;
+    // Loop over tracks
+    auto chargedTracks = chargedTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
+    float ux = 0., uy = 0., v1a = 0., v1c = 0.;
+    for (auto const& track : chargedTracks) {
+      // Fill track QA
+      histos.fill(HIST("TrackQA/hPtDcaZ"), track.pt(), track.dcaZ());
+      histos.fill(HIST("TrackQA/hPtDcaXY"), track.pt(), track.dcaXY());
+      histos.fill(HIST("TrackQA/hTrackTPCdEdX"), track.pt(), track.tpcSignal());
+
+      // Get directed flow
+      ux = std::cos(track.phi());
+      uy = std::sin(track.phi());
+      v1a = ux * vSP[kXa] + uy * vSP[kYa];
+      v1c = ux * vSP[kXc] + uy * vSP[kYc];
+
+      // Charged particle directed flow
+      histos.fill(HIST("DF/hAQu"), cent, track.eta(), v1a);
+      histos.fill(HIST("DF/hCQu"), cent, track.eta(), v1c);
+      if (track.sign() > 0) {
+        histos.fill(HIST("DF/hAQuPos"), cent, track.eta(), v1a);
+        histos.fill(HIST("DF/hCQuPos"), cent, track.eta(), v1c);
+      } else {
+        histos.fill(HIST("DF/hAQuNeg"), cent, track.eta(), v1a);
+        histos.fill(HIST("DF/hCQuNeg"), cent, track.eta(), v1c);
+      }
     }
+  }
 
-    auto zdc = bc.zdc();
-    auto znaEnergy = zdc.energySectorZNA();
-    auto zncEnergy = zdc.energySectorZNC();
-    auto znaEnergyCommon = zdc.energyCommonZNA();
-    auto zncEnergyCommon = zdc.energyCommonZNC();
+  PROCESS_SWITCH(FlowEventPlane, processChargedFlow, "Charged particle flow process", false);
 
-    // check energy deposits
-    if (znaEnergyCommon <= 0 || zncEnergyCommon <= 0 || znaEnergy[0] <= 0 || znaEnergy[1] <= 0 || znaEnergy[2] <= 0 || znaEnergy[3] <= 0 || zncEnergy[0] <= 0 || zncEnergy[1] <= 0 || zncEnergy[2] <= 0 || zncEnergy[3] <= 0) {
+  void processIdFlow(CollisionsRun3::iterator const& collision, Tracks const&)
+  {
+    if (!selCollision(collision, vSP)) {
       return;
     }
+    // Loop over tracks
+    auto pionTracks = pionTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
+    auto kaonTracks = kaonTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
+    auto protonTracks = protonTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
+    getIdFlow<kPi>(pionTracks);
+    getIdFlow<kKa>(kaonTracks);
+    getIdFlow<kPr>(protonTracks);
+  }
 
-    // Fill collision table
-    histos.fill(HIST("hCent"), cent);
-    colSpTable(runNum, timestamp, cent, posX, posY, posZ, znaEnergyCommon, zncEnergyCommon, znaEnergy.data(), zncEnergy.data());
-    colIdx = colSpTable.lastIndex();
+  PROCESS_SWITCH(FlowEventPlane, processIdFlow, "Identified particle flow process", false);
 
-    // Analyze Tracks [charged, pi, K, p]
-    analyzeIdHadrons(tracks);
+  void processResoFlow(CollisionsRun3::iterator const& collision, Tracks const&)
+  {
+    if (!selCollision(collision, vSP)) {
+      return;
+    }
 
-    // Analyze V0s [Lambda, K0S]
-    analyzeV0s(V0s, tracks);
+    // Track partitions
+    auto kaonTracks = kaonTrackPartition->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
 
-    // Done
-    return;
+    // Resonance flow
+    getResoFlow<kPhi0>(kaonTracks, kaonTracks, vSP);
   }
+  PROCESS_SWITCH(FlowEventPlane, processResoFlow, "Resonance flow process", false);
 
-  using BCsRun3 = soa::Join<aod::BCsWithTimestamps, aod::Run3MatchedToBCSparse>;
-  using CollisionsRun3 = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::CentFT0Ms, aod::CentFV0As, aod::TPCMults, aod::PVMults, aod::MultsGlobal, aod::MultsExtra>;
-  using Tracks = soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::TOFSignal, aod::pidTOFbeta, aod::pidTPCEl, aod::pidTPCPi, aod::pidTOFPi, aod::pidTPCKa, aod::pidTOFKa, aod::pidTPCPr, aod::pidTOFPr, aod::TrackCompColls>;
+  void processV0Flow(CollisionsRun3::iterator const& collision, aod::V0Datas const& V0s, TracksV0s const& tracks)
+  {
+    if (!selCollision(collision, vSP)) {
+      return;
+    }
 
-  void processDummy(CollisionsRun3::iterator const&) {}
+    // Loop over v0s
+    for (auto const& v0 : V0s) {
+      // Topological and kinematic selections
+      if (std::abs(v0.eta()) >= cEtaCut || v0.dcaV0daughters() >= cDcaV0Dau || v0.v0radius() <= cMinV0Radius || v0.v0Type() != cV0TypeSelection) {
+        continue;
+      }
 
-  PROCESS_SWITCH(FlowEventPlane, processDummy, "Dummy process", true);
+      // Directed flow
+      float ux = std::cos(v0.phi());
+      float uy = std::sin(v0.phi());
+      float v1a = ux * vSP[kXa] + uy * vSP[kYa];
+      float v1c = ux * vSP[kXc] + uy * vSP[kYc];
+
+      // Initialize selection objects
+      auto postrack = v0.template posTrack_as<TracksV0s>();
+      auto negtrack = v0.template negTrack_as<TracksV0s>();
+      float ctauK0Short = v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * MassKaonNeutral;
+      float ctauLambda = v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * MassLambda0;
+
+      // K0Short
+      if (selV0DauTracks<kK0S>(v0, postrack, negtrack) && v0.v0cosPA() > cK0ShortCosPA && ctauK0Short < cK0ShortCTau && v0.qtarm() >= cArmPodSel * std::abs(v0.alpha()) && v0.pt() >= cK0SMinPt && v0.pt() < cK0SMaxPt) {
+        fillV0QAHist<kK0S>(collision, v0, tracks);
+        histos.fill(HIST("V0/K0Short/hMassVsRap"), cent, v0.eta(), v0.mK0Short());
+        histos.fill(HIST("V0/K0Short/Flow/hQuA"), cent, v0.eta(), v0.mK0Short(), v1a);
+        histos.fill(HIST("V0/K0Short/Flow/hQuC"), cent, v0.eta(), v0.mK0Short(), v1c);
+      }
 
-  void processFEP(CollisionsRun3::iterator const& collision, BCsRun3 const&, aod::Zdcs const&,
-                  Tracks const& tracks, aod::V0Datas const& v0s)
-  {
-    // Get bunch crossing
-    auto bc = collision.template foundBC_as<BCsRun3>();
-    runNum = collision.template foundBC_as<BCsRun3>().runNumber();
-    timestamp = collision.template foundBC_as<BCsRun3>().timestamp();
-    analyzeCollision(bc, collision, tracks, v0s);
-  }
+      // Lambda
+      if (selV0DauTracks<kLambda>(v0, postrack, negtrack) && v0.v0cosPA() > cLambdaCosPA && ctauLambda < cLambdaCTau && std::abs(v0.mK0Short() - MassK0Short) >= cK0SMassRej && v0.qtarm() < cArmPodSel * std::abs(v0.alpha()) && v0.alpha() > 0 && v0.pt() >= cLambdaMinPt && v0.pt() < cLambdaMaxPt) {
+        fillV0QAHist<kLambda>(collision, v0, tracks);
+        histos.fill(HIST("V0/Lambda/hMassVsRap"), cent, v0.eta(), v0.mLambda());
+        histos.fill(HIST("V0/Lambda/Flow/hQuA"), cent, v0.eta(), v0.mLambda(), v1a);
+        histos.fill(HIST("V0/Lambda/Flow/hQuC"), cent, v0.eta(), v0.mLambda(), v1c);
+        histos.fill(HIST("V0/LambdaAntiLambda/hMassVsRap"), cent, v0.eta(), v0.mLambda());
+        histos.fill(HIST("V0/LambdaAntiLambda/Flow/hQuA"), cent, v0.eta(), v0.mLambda(), v1a);
+        histos.fill(HIST("V0/LambdaAntiLambda/Flow/hQuC"), cent, v0.eta(), v0.mLambda(), v1c);
+      }
 
-  PROCESS_SWITCH(FlowEventPlane, processFEP, "Flow Event Plane Table Producer", false);
+      // AntiLambda
+      if (selV0DauTracks<kAntiLambda>(v0, postrack, negtrack) && v0.v0cosPA() > cLambdaCosPA && ctauLambda < cLambdaCTau && std::abs(v0.mK0Short() - MassK0Short) >= cK0SMassRej && v0.qtarm() < cArmPodSel * std::abs(v0.alpha()) && v0.alpha() < 0 && v0.pt() >= cLambdaMinPt && v0.pt() < cLambdaMaxPt) {
+        fillV0QAHist<kAntiLambda>(collision, v0, tracks);
+        histos.fill(HIST("V0/AntiLambda/hMassVsRap"), cent, v0.eta(), v0.mAntiLambda());
+        histos.fill(HIST("V0/AntiLambda/Flow/hQuA"), cent, v0.eta(), v0.mAntiLambda(), v1a);
+        histos.fill(HIST("V0/AntiLambda/Flow/hQuC"), cent, v0.eta(), v0.mAntiLambda(), v1c);
+        histos.fill(HIST("V0/LambdaAntiLambda/hMassVsRap"), cent, v0.eta(), v0.mAntiLambda());
+        histos.fill(HIST("V0/LambdaAntiLambda/Flow/hQuA"), cent, v0.eta(), v0.mAntiLambda(), v1a);
+        histos.fill(HIST("V0/LambdaAntiLambda/Flow/hQuC"), cent, v0.eta(), v0.mAntiLambda(), v1c);
+      }
+    }
+  }
+  PROCESS_SWITCH(FlowEventPlane, processV0Flow, "Lambda flow process", false);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
   return WorkflowSpec{
+    adaptAnalysisTask<SpCalibTableProducer>(cfgc),
+    adaptAnalysisTask<SpectatorPlaneTableProducer>(cfgc),
     adaptAnalysisTask<FlowEventPlane>(cfgc)};
 }