diff --git a/PWGLF/TableProducer/Strangeness/cascadeflow.cxx b/PWGLF/TableProducer/Strangeness/cascadeflow.cxx
index c66a1db16fd..5179a89846b 100644
--- a/PWGLF/TableProducer/Strangeness/cascadeflow.cxx
+++ b/PWGLF/TableProducer/Strangeness/cascadeflow.cxx
@@ -9,8 +9,11 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 ///
+/// \file cascadeflow.cxx
+///
 /// \brief Task to create derived data for cascade flow analyses
-/// \authors: Chiara De Martin (chiara.de.martin@cern.ch), Maximiliano Puccio (maximiliano.puccio@cern.ch)
+/// \author Chiara De Martin (chiara.de.martin@cern.ch)
+/// \author Maximiliano Puccio (maximiliano.puccio@cern.ch)
 
 #include <vector>
 #include <string>
@@ -40,17 +43,21 @@ using std::array;
 
 using DauTracks = soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>;
 using CollEventPlane = soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraFT0CQVs, aod::StraFT0CQVsEv, aod::StraTPCQVs>::iterator;
+using CollEventPlaneCentralFW = soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraFT0CQVs, aod::StraTPCQVs>::iterator;
 using CollEventAndSpecPlane = soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraFT0CQVs, aod::StraFT0CQVsEv, aod::StraTPCQVs, aod::StraZDCSP>::iterator;
-using CollEventPlaneCentralFW = soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraFT0CQVs, aod::StraTPCQVs, aod::StraZDCSP>::iterator;
+using CollEventAndSpecPlaneCentralFW = soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraFT0CQVs, aod::StraTPCQVs, aod::StraZDCSP>::iterator;
 using MCCollisionsStra = soa::Join<aod::StraMCCollisions, aod::StraMCCollMults>;
 using CascCandidates = soa::Join<aod::CascCollRefs, aod::CascCores, aod::CascExtras, aod::CascBBs>;
 using CascMCCandidates = soa::Join<aod::CascCollRefs, aod::CascCores, aod::CascExtras, aod::CascBBs, aod::CascCoreMCLabels>;
 
+const int nParticles = 2;
+const int nParameters = 4;
+
 namespace cascadev2
 {
 enum species { Xi = 0,
                Omega = 1 };
-constexpr double massSigmaParameters[4][2]{
+constexpr double massSigmaParameters[nParameters][nParticles]{
   {4.9736e-3, 0.006815},
   {-2.39594, -2.257},
   {1.8064e-3, 0.00138},
@@ -61,13 +68,13 @@ const double AlphaXi[2] = {-0.390, 0.371};     // decay parameter of XiMinus and
 const double AlphaOmega[2] = {0.0154, -0.018}; // decay parameter of OmegaMinus and OmegaPlus
 const double AlphaLambda[2] = {0.747, -0.757}; // decay parameter of Lambda and AntiLambda
 
-std::shared_ptr<TH2> hMassBeforeSelVsPt[2];
-std::shared_ptr<TH2> hMassAfterSelVsPt[2];
-std::shared_ptr<TH1> hSignalScoreBeforeSel[2];
-std::shared_ptr<TH1> hBkgScoreBeforeSel[2];
-std::shared_ptr<TH1> hSignalScoreAfterSel[2];
-std::shared_ptr<TH1> hBkgScoreAfterSel[2];
-std::shared_ptr<THn> hSparseV2C[2];
+std::shared_ptr<TH2> hMassBeforeSelVsPt[nParticles];
+std::shared_ptr<TH2> hMassAfterSelVsPt[nParticles];
+std::shared_ptr<TH1> hSignalScoreBeforeSel[nParticles];
+std::shared_ptr<TH1> hBkgScoreBeforeSel[nParticles];
+std::shared_ptr<TH1> hSignalScoreAfterSel[nParticles];
+std::shared_ptr<TH1> hBkgScoreAfterSel[nParticles];
+std::shared_ptr<THn> hSparseV2C[nParticles];
 } // namespace cascadev2
 
 namespace cascade_flow_cuts_ml
@@ -154,7 +161,7 @@ struct cascadeFlow {
   ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {VARIABLE_WIDTH, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5, 6, 8, 10}, "#it{p}_{T} (GeV/#it{c})"};
   ConfigurableAxis thnConfigAxisPtLambda{"thnConfigAxisPtLambda", {VARIABLE_WIDTH, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5, 6, 8, 10, 20}, "#it{p}_{T} (GeV/#it{c})"};
   ConfigurableAxis thnConfigAxisCharge{"thnConfigAxisCharge", {2, 0, 2}, ""};
-  ConfigurableAxis thnConfigAxisPsiDiff{"thnConfigAxisPsiDiff", {100, 0, 2 * TMath::Pi()}, ""};
+  ConfigurableAxis thnConfigAxisPsiDiff{"thnConfigAxisPsiDiff", {100, 0, o2::constants::math::TwoPI}, ""};
   ConfigurableAxis thnConfigAxisMassXi{"thnConfigAxisMassXi", {45, 1.300, 1.345}, ""};
   ConfigurableAxis thnConfigAxisMassOmega{"thnConfigAxisMassOmega", {45, 1.655, 1.690}, ""};
   ConfigurableAxis thnConfigAxisMassLambda{"thnConfigAxisMassLambda", {60, 1.1, 1.13}, ""};
@@ -261,7 +268,7 @@ struct cascadeFlow {
       histos.fill(HIST("hNEvents"), 1.5);
 
     // Z vertex selection
-    if (TMath::Abs(collision.posZ()) > cutzvertex) {
+    if (std::abs(collision.posZ()) > cutzvertex) {
       return false;
     }
     if (isFillHisto)
@@ -356,10 +363,10 @@ struct cascadeFlow {
   double GetPhiInRange(double phi)
   {
     while (phi < 0) {
-      phi += TMath::Pi();
+      phi += o2::constants::math::PI;
     }
-    while (phi > TMath::Pi()) {
-      phi -= TMath::Pi();
+    while (phi > o2::constants::math::PI) {
+      phi -= o2::constants::math::PI;
     }
     return phi;
   }
@@ -378,13 +385,13 @@ struct cascadeFlow {
   Produces<aod::CascAnalysis> analysisSample;
   Configurable<LabeledArray<double>> parSigmaMass{
     "parSigmaMass",
-    {cascadev2::massSigmaParameters[0], 4, 2,
+    {cascadev2::massSigmaParameters[0], nParameters, nParticles,
      cascadev2::massSigmaParameterNames, cascadev2::speciesNames},
     "Mass resolution parameters: [0]*exp([1]*x)+[2]*exp([3]*x)"};
 
   float getNsigmaMass(const cascadev2::species s, const float pt, const float nsigma = 6)
   {
-    const auto sigma = parSigmaMass->get(0u, s) * exp(parSigmaMass->get(1, s) * pt) + parSigmaMass->get(2, s) * exp(parSigmaMass->get(3, s) * pt);
+    const auto sigma = parSigmaMass->get(0u, s) * std::exp(parSigmaMass->get(1, s) * pt) + parSigmaMass->get(2, s) * std::exp(parSigmaMass->get(3, s) * pt);
     return nsigma * sigma;
   }
 
@@ -416,9 +423,9 @@ struct cascadeFlow {
   template <class collision_t, class cascade_t>
   void fillAnalysedTable(collision_t coll, bool hasEventPlane, bool hasSpectatorPlane, cascade_t casc, float v2CSP, float v2CEP, float v1SP_ZDCA, float v1SP_ZDCC, float PsiT0C, float BDTresponseXi, float BDTresponseOmega, int pdgCode)
   {
-    double masses[2]{o2::constants::physics::MassXiMinus, o2::constants::physics::MassOmegaMinus};
-    ROOT::Math::PxPyPzMVector cascadeVector[2], lambdaVector, protonVector;
-    float cosThetaStarLambda[2], cosThetaStarProton;
+    double masses[nParticles]{o2::constants::physics::MassXiMinus, o2::constants::physics::MassOmegaMinus};
+    ROOT::Math::PxPyPzMVector cascadeVector[nParticles], lambdaVector, protonVector;
+    float cosThetaStarLambda[nParticles], cosThetaStarProton;
     lambdaVector.SetCoordinates(casc.pxlambda(), casc.pylambda(), casc.pzlambda(), o2::constants::physics::MassLambda);
     ROOT::Math::Boost lambdaBoost{lambdaVector.BoostToCM()};
     if (casc.sign() > 0) {
@@ -428,7 +435,7 @@ struct cascadeFlow {
     }
     auto boostedProton{lambdaBoost(protonVector)};
     cosThetaStarProton = boostedProton.Pz() / boostedProton.P();
-    for (int i{0}; i < 2; ++i) {
+    for (int i{0}; i < nParticles; ++i) {
       cascadeVector[i].SetCoordinates(casc.px(), casc.py(), casc.pz(), masses[i]);
       ROOT::Math::Boost cascadeBoost{cascadeVector[i].BoostToCM()};
       auto boostedLambda{cascadeBoost(lambdaVector)};
@@ -531,9 +538,9 @@ struct cascadeFlow {
     }
     histos.add("hEventVertexZ", "hEventVertexZ", kTH1F, {{120, -12., 12.}});
     histos.add("hEventCentrality", "hEventCentrality", kTH1F, {{101, 0, 101}});
-    histos.add("hPsiT0C", "hPsiT0C", HistType::kTH1D, {{100, -TMath::Pi(), TMath::Pi()}});
-    histos.add("hPsiT0CvsCentFT0C", "hPsiT0CvsCentFT0C", HistType::kTH2D, {CentAxis, {100, -TMath::Pi(), TMath::Pi()}});
-    histos.add("hPsiZDCA_vs_ZDCC", "hPsiZDCA_vs_ZDCC", HistType::kTH2D, {{100, -TMath::Pi(), TMath::Pi()}, {100, -TMath::Pi(), TMath::Pi()}});
+    histos.add("hPsiT0C", "hPsiT0C", HistType::kTH1D, {{100, -o2::constants::math::PI, o2::constants::math::PI}});
+    histos.add("hPsiT0CvsCentFT0C", "hPsiT0CvsCentFT0C", HistType::kTH2D, {CentAxis, {100, -o2::constants::math::PI, o2::constants::math::PI}});
+    histos.add("hPsiZDCA_vs_ZDCC", "hPsiZDCA_vs_ZDCC", HistType::kTH2D, {{100, -o2::constants::math::PI, o2::constants::math::PI}, {100, -o2::constants::math::PI, o2::constants::math::PI}});
     histos.add("hEventNchCorrelation", "hEventNchCorrelation", kTH2F, {{5000, 0, 5000}, {5000, 0, 2500}});
     histos.add("hEventPVcontributorsVsCentrality", "hEventPVcontributorsVsCentrality", kTH2F, {{100, 0, 100}, {5000, 0, 5000}});
     histos.add("hEventGlobalTracksVsCentrality", "hEventGlobalTracksVsCentrality", kTH2F, {{100, 0, 100}, {2500, 0, 2500}});
@@ -554,8 +561,8 @@ struct cascadeFlow {
     histos.add("hOmegaPtvsCent", "hOmegaPtvsCent", HistType::kTH2F, {{100, 0, 100}, {400, 0, 20}});
     histos.add("hOmegaPtvsCentEta08", "hOmegaPtvsCentEta08", HistType::kTH2F, {{100, 0, 100}, {400, 0, 20}});
     histos.add("hOmegaPtvsCentY05", "hOmegaPtvsCentY05", HistType::kTH2F, {{100, 0, 100}, {400, 0, 20}});
-    histos.add("hCascadePhi", "hCascadePhi", HistType::kTH1F, {{100, 0, 2 * TMath::Pi()}});
-    histos.add("hcascminuspsiT0C", "hcascminuspsiT0C", HistType::kTH1F, {{100, 0, TMath::Pi()}});
+    histos.add("hCascadePhi", "hCascadePhi", HistType::kTH1F, {{100, 0, o2::constants::math::TwoPI}});
+    histos.add("hcascminuspsiT0C", "hcascminuspsiT0C", HistType::kTH1F, {{100, 0, o2::constants::math::PI}});
     histos.add("hv2CEPvsFT0C", "hv2CEPvsFT0C", HistType::kTH2F, {CentAxis, {100, -1, 1}});
     histos.add("hv2CEPvsv2CSP", "hv2CEPvsV2CSP", HistType::kTH2F, {{100, -1, 1}, {100, -1, 1}});
     histos.add("hv1EPvsv1SP", "hV1EPvsV1SP", HistType::kTH2F, {{100, -1, 1}, {100, -1, 1}});
@@ -653,7 +660,7 @@ struct cascadeFlow {
       histosMCGen.get<TH1>(HIST("hNCascGen"))->GetXaxis()->SetBinLabel(n, hNCascLabelsMC[n - 1]);
     }
 
-    for (int iS{0}; iS < 2; ++iS) {
+    for (int iS{0}; iS < nParticles; ++iS) {
       cascadev2::hMassBeforeSelVsPt[iS] = histos.add<TH2>(Form("hMassBeforeSelVsPt%s", cascadev2::speciesNames[iS].data()), "hMassBeforeSelVsPt", HistType::kTH2F, {massCascAxis[iS], ptAxis});
       cascadev2::hMassAfterSelVsPt[iS] = histos.add<TH2>(Form("hMassAfterSelVsPt%s", cascadev2::speciesNames[iS].data()), "hMassAfterSelVsPt", HistType::kTH2F, {massCascAxis[iS], ptAxis});
       cascadev2::hSignalScoreBeforeSel[iS] = histos.add<TH1>(Form("hSignalScoreBeforeSel%s", cascadev2::speciesNames[iS].data()), "Signal score before selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}});
@@ -699,7 +706,7 @@ struct cascadeFlow {
     histos.fill(HIST("hEventCentrality"), coll.centFT0C());
     histos.fill(HIST("hEventVertexZ"), coll.posZ());
 
-    for (auto& casc : Cascades) {
+    for (auto const& casc : Cascades) {
       if (gRandom->Uniform() > downsample) {
         continue;
       }
@@ -749,14 +756,14 @@ struct cascadeFlow {
     histos.fill(HIST("hEventCentrality"), coll.centFT0C());
     histos.fill(HIST("hEventVertexZ"), coll.posZ());
 
-    for (auto& casc : Cascades) {
+    for (auto const& casc : Cascades) {
       if (!casc.has_cascMCCore())
         continue;
 
       auto cascMC = casc.cascMCCore_as<soa::Join<aod::CascMCCores, aod::CascMCCollRefs>>();
       int pdgCode{cascMC.pdgCode()};
-      if (!(std::abs(pdgCode) == 3312 && std::abs(cascMC.pdgCodeV0()) == 3122 && std::abs(cascMC.pdgCodeBachelor()) == 211)     // Xi
-          && !(std::abs(pdgCode) == 3334 && std::abs(cascMC.pdgCodeV0()) == 3122 && std::abs(cascMC.pdgCodeBachelor()) == 321)) // Omega
+      if (!(std::abs(pdgCode) == PDG_t::kXiMinus && std::abs(cascMC.pdgCodeV0()) == PDG_t::kLambda0 && std::abs(cascMC.pdgCodeBachelor()) == PDG_t::kPiPlus) // Xi
+          && !(std::abs(pdgCode) == PDG_t::kOmegaMinus && std::abs(cascMC.pdgCodeV0()) == PDG_t::kLambda0 && std::abs(cascMC.pdgCodeBachelor()) == kKPlus))  // Omega
         continue;
 
       auto negExtra = casc.negTrackExtra_as<DauTracks>();
@@ -821,8 +828,8 @@ struct cascadeFlow {
     resolution.fill(HIST("QVectorsNormTPCAC"), eventplaneVecTPCA.Dot(eventplaneVecTPCC) / (coll.qTPCR() * coll.qTPCL()), coll.centFT0C());
     resolution.fill(HIST("QVectorsSpecPlane"), spectatorplaneVecZDCC.Dot(spectatorplaneVecZDCA), coll.centFT0C());
 
-    std::vector<float> bdtScore[2];
-    for (auto& casc : Cascades) {
+    std::vector<float> bdtScore[nParticles];
+    for (auto const& casc : Cascades) {
 
       /// Add some minimal cuts for single track variables (min number of TPC clusters)
       auto negExtra = casc.negTrackExtra_as<DauTracks>();
@@ -864,7 +871,7 @@ struct cascadeFlow {
         isSelectedCasc[0] = mlResponseXi.isSelectedMl(inputFeaturesCasc, casc.pt(), bdtScore[0]);
         isSelectedCasc[1] = mlResponseOmega.isSelectedMl(inputFeaturesCasc, casc.pt(), bdtScore[1]);
 
-        for (int iS{0}; iS < 2; ++iS) {
+        for (int iS{0}; iS < nParticles; ++iS) {
           // Fill BDT score histograms before selection
           cascadev2::hSignalScoreBeforeSel[iS]->Fill(bdtScore[0][1]);
           cascadev2::hBkgScoreBeforeSel[iS]->Fill(bdtScore[1][0]);
@@ -884,12 +891,12 @@ struct cascadeFlow {
       ROOT::Math::XYZVector cascQvec{std::cos(2 * casc.phi()), std::sin(2 * casc.phi()), 0};
       auto v2CSP = cascQvec.Dot(eventplaneVecT0C); // not normalised by amplitude
       auto cascminuspsiT0C = GetPhiInRange(casc.phi() - PsiT0C);
-      auto v2CEP = TMath::Cos(2.0 * cascminuspsiT0C);
+      auto v2CEP = std::cos(2.0 * cascminuspsiT0C);
       ROOT::Math::XYZVector cascUvec{std::cos(casc.phi()), std::sin(casc.phi()), 0};
       auto v1SP_ZDCA = cascUvec.Dot(spectatorplaneVecZDCA);
       auto v1SP_ZDCC = cascUvec.Dot(spectatorplaneVecZDCC);
-      auto v1EP_ZDCA = TMath::Cos(casc.phi() - coll.psiZDCA());
-      auto v1EP_ZDCC = TMath::Cos(casc.phi() - coll.psiZDCC());
+      auto v1EP_ZDCA = std::cos(casc.phi() - coll.psiZDCA());
+      auto v1EP_ZDCC = std::cos(casc.phi() - coll.psiZDCC());
       float v1SP = 0.5 * (v1SP_ZDCA - v1SP_ZDCC);
       float v1EP = 0.5 * (v1EP_ZDCA - v1EP_ZDCC); // same as v1SP
 
@@ -906,14 +913,14 @@ struct cascadeFlow {
       }
       auto boostedProton{lambdaBoost(protonVector)};
       cosThetaStarProton = boostedProton.Pz() / boostedProton.P();
-      for (int i{0}; i < 2; ++i) {
+      for (int i{0}; i < nParticles; ++i) {
         cascadeVector[i].SetCoordinates(casc.px(), casc.py(), casc.pz(), masses[i]);
         ROOT::Math::Boost cascadeBoost{cascadeVector[i].BoostToCM()};
         auto boostedLambda{cascadeBoost(lambdaVector)};
         cosThetaStarLambda[i] = boostedLambda.Pz() / boostedLambda.P();
       }
 
-      double ptLambda = sqrt(pow(casc.pxlambda(), 2) + pow(casc.pylambda(), 2));
+      double ptLambda = std::sqrt(std::pow(casc.pxlambda(), 2) + std::pow(casc.pylambda(), 2));
       auto etaLambda = RecoDecay::eta(std::array{casc.pxlambda(), casc.pylambda(), casc.pzlambda()});
 
       // acceptance values if requested
@@ -966,7 +973,7 @@ struct cascadeFlow {
         cascadev2::hSparseV2C[0]->Fill(values);
       }
 
-      float BDTresponse[2]{0.f, 0.f};
+      float BDTresponse[nParticles]{0.f, 0.f};
       if (isApplyML) {
         BDTresponse[0] = bdtScore[0][1];
         BDTresponse[1] = bdtScore[1][1];
@@ -1054,7 +1061,7 @@ struct cascadeFlow {
     }
   }
 
-  void processAnalyseDataEPCentralFW(CollEventPlaneCentralFW const& coll, CascCandidates const& Cascades, DauTracks const&)
+  void processAnalyseDataEP2CentralFW(CollEventPlaneCentralFW const& coll, CascCandidates const& Cascades, DauTracks const&)
   {
 
     if (!AcceptEvent(coll, 1)) {
@@ -1063,14 +1070,281 @@ struct cascadeFlow {
 
     // select only events used for the calibration of the event plane
     if (isGoodEventEP) {
-      if (abs(coll.qvecFT0CRe()) > 990 || abs(coll.qvecFT0CIm()) > 990 || abs(coll.qvecBNegRe()) > 990 || abs(coll.qvecBNegIm()) > 990 || abs(coll.qvecBPosRe()) > 990 || abs(coll.qvecBPosIm()) > 990) {
+      if (std::abs(coll.qvecFT0CRe()) > 990 || std::abs(coll.qvecFT0CIm()) > 990 || std::abs(coll.qvecBNegRe()) > 990 || std::abs(coll.qvecBNegIm()) > 990 || std::abs(coll.qvecBPosRe()) > 990 || std::abs(coll.qvecBPosIm()) > 990) {
         return;
       }
     }
 
-    // event has event plane
+    // event has FT0C event plane
+    bool hasEventPlane = 0;
+    if (std::abs(coll.qvecFT0CRe()) < 990 && std::abs(coll.qvecFT0CIm()) < 990)
+      hasEventPlane = 1;
+
+    histos.fill(HIST("hNEvents"), 9.5);
+    histos.fill(HIST("hEventNchCorrelationAfterEP"), coll.multNTracksPVeta1(), coll.multNTracksGlobal());
+    histos.fill(HIST("hEventPVcontributorsVsCentralityAfterEP"), coll.centFT0C(), coll.multNTracksPVeta1());
+    histos.fill(HIST("hEventGlobalTracksVsCentralityAfterEP"), coll.centFT0C(), coll.multNTracksGlobal());
+
+    histos.fill(HIST("hEventCentrality"), coll.centFT0C());
+    histos.fill(HIST("hEventVertexZ"), coll.posZ());
+
+    ROOT::Math::XYZVector eventplaneVecT0C{coll.qvecFT0CRe(), coll.qvecFT0CIm(), 0};
+    ROOT::Math::XYZVector eventplaneVecTPCA{coll.qvecBPosRe(), coll.qvecBPosIm(), 0};
+    ROOT::Math::XYZVector eventplaneVecTPCC{coll.qvecBNegRe(), coll.qvecBNegIm(), 0};
+
+    const float PsiT0C = std::atan2(coll.qvecFT0CIm(), coll.qvecFT0CRe()) * 0.5f;
+    histos.fill(HIST("hPsiT0C"), PsiT0C);
+    histos.fill(HIST("hPsiT0CvsCentFT0C"), coll.centFT0C(), PsiT0C);
+
+    resolution.fill(HIST("QVectorsT0CTPCA"), eventplaneVecT0C.Dot(eventplaneVecTPCA), coll.centFT0C());
+    resolution.fill(HIST("QVectorsT0CTPCC"), eventplaneVecT0C.Dot(eventplaneVecTPCC), coll.centFT0C());
+    resolution.fill(HIST("QVectorsTPCAC"), eventplaneVecTPCA.Dot(eventplaneVecTPCC), coll.centFT0C());
+    resolution.fill(HIST("QVectorsNormT0CTPCA"), eventplaneVecT0C.Dot(eventplaneVecTPCA) / (coll.qTPCR() * coll.sumAmplFT0C()), coll.centFT0C());
+    resolution.fill(HIST("QVectorsNormT0CTPCC"), eventplaneVecT0C.Dot(eventplaneVecTPCC) / (coll.qTPCL() * coll.sumAmplFT0C()), coll.centFT0C());
+    resolution.fill(HIST("QVectorsNormTPCAC"), eventplaneVecTPCA.Dot(eventplaneVecTPCC) / (coll.qTPCR() * coll.qTPCL()), coll.centFT0C());
+
+    std::vector<float> bdtScore[nParticles];
+    for (auto const& casc : Cascades) {
+
+      /// Add some minimal cuts for single track variables (min number of TPC clusters)
+      auto negExtra = casc.negTrackExtra_as<DauTracks>();
+      auto posExtra = casc.posTrackExtra_as<DauTracks>();
+      auto bachExtra = casc.bachTrackExtra_as<DauTracks>();
+
+      int counter = 0;
+      IsCascAccepted(casc, negExtra, posExtra, bachExtra, counter);
+      histos.fill(HIST("hCascade"), counter);
+
+      // ML selections
+      bool isSelectedCasc[nParticles]{false, false};
+
+      std::vector<float> inputFeaturesCasc{casc.cascradius(),
+                                           casc.v0radius(),
+                                           casc.casccosPA(coll.posX(), coll.posY(), coll.posZ()),
+                                           casc.v0cosPA(coll.posX(), coll.posY(), coll.posZ()),
+                                           casc.dcapostopv(),
+                                           casc.dcanegtopv(),
+                                           casc.dcabachtopv(),
+                                           casc.dcacascdaughters(),
+                                           casc.dcaV0daughters(),
+                                           casc.dcav0topv(coll.posX(), coll.posY(), coll.posZ()),
+                                           casc.bachBaryonCosPA(),
+                                           casc.bachBaryonDCAxyToPV()};
+
+      float massCasc[nParticles]{casc.mXi(), casc.mOmega()};
+
+      // pt cut
+      if (casc.pt() < CandidateConfigs.MinPt || casc.pt() > CandidateConfigs.MaxPt) {
+        continue;
+      }
+
+      cascadev2::hMassBeforeSelVsPt[0]->Fill(massCasc[0], casc.pt());
+      cascadev2::hMassBeforeSelVsPt[1]->Fill(massCasc[1], casc.pt());
+
+      if (isApplyML) {
+        // Retrieve model output and selection outcome
+        isSelectedCasc[0] = mlResponseXi.isSelectedMl(inputFeaturesCasc, casc.pt(), bdtScore[0]);
+        isSelectedCasc[1] = mlResponseOmega.isSelectedMl(inputFeaturesCasc, casc.pt(), bdtScore[1]);
+
+        for (int iS{0}; iS < nParticles; ++iS) {
+          // Fill BDT score histograms before selection
+          cascadev2::hSignalScoreBeforeSel[iS]->Fill(bdtScore[0][1]);
+          cascadev2::hBkgScoreBeforeSel[iS]->Fill(bdtScore[1][0]);
+
+          // Fill histograms for selected candidates
+          if (isSelectedCasc[iS]) {
+            cascadev2::hSignalScoreAfterSel[iS]->Fill(bdtScore[0][1]);
+            cascadev2::hBkgScoreAfterSel[iS]->Fill(bdtScore[1][0]);
+            cascadev2::hMassAfterSelVsPt[iS]->Fill(massCasc[iS], casc.pt());
+          }
+        }
+      } else {
+        isSelectedCasc[0] = true;
+        isSelectedCasc[1] = true;
+      }
+
+      ROOT::Math::XYZVector cascQvec{std::cos(2 * casc.phi()), std::sin(2 * casc.phi()), 0};
+      auto v2CSP = cascQvec.Dot(eventplaneVecT0C); // not normalised by amplitude
+      auto cascminuspsiT0C = GetPhiInRange(casc.phi() - PsiT0C);
+      auto v2CEP = std::cos(2.0 * cascminuspsiT0C);
+      ROOT::Math::XYZVector cascUvec{std::cos(casc.phi()), std::sin(casc.phi()), 0};
+
+      // polarization variables
+      double masses[nParticles]{o2::constants::physics::MassXiMinus, o2::constants::physics::MassOmegaMinus};
+      ROOT::Math::PxPyPzMVector cascadeVector[nParticles], lambdaVector, protonVector;
+      float cosThetaStarLambda[nParticles], cosThetaStarProton;
+      lambdaVector.SetCoordinates(casc.pxlambda(), casc.pylambda(), casc.pzlambda(), o2::constants::physics::MassLambda);
+      ROOT::Math::Boost lambdaBoost{lambdaVector.BoostToCM()};
+      if (casc.sign() > 0) {
+        protonVector.SetCoordinates(casc.pxneg(), casc.pyneg(), casc.pzneg(), o2::constants::physics::MassProton);
+      } else {
+        protonVector.SetCoordinates(casc.pxpos(), casc.pypos(), casc.pzpos(), o2::constants::physics::MassProton);
+      }
+      auto boostedProton{lambdaBoost(protonVector)};
+      cosThetaStarProton = boostedProton.Pz() / boostedProton.P();
+      for (int i{0}; i < nParticles; ++i) {
+        cascadeVector[i].SetCoordinates(casc.px(), casc.py(), casc.pz(), masses[i]);
+        ROOT::Math::Boost cascadeBoost{cascadeVector[i].BoostToCM()};
+        auto boostedLambda{cascadeBoost(lambdaVector)};
+        cosThetaStarLambda[i] = boostedLambda.Pz() / boostedLambda.P();
+      }
+
+      double ptLambda = std::sqrt(std::pow(casc.pxlambda(), 2) + std::pow(casc.pylambda(), 2));
+      auto etaLambda = RecoDecay::eta(std::array{casc.pxlambda(), casc.pylambda(), casc.pzlambda()});
+
+      // acceptance values if requested
+      double MeanCos2ThetaLambdaFromXi = 1;
+      double MeanCos2ThetaLambdaFromOmega = 1;
+      double MeanCos2ThetaProtonFromLambda = 1;
+      if (applyAcceptanceCorrection) {
+        if (ptLambda < CandidateConfigs.MinPtLambda || ptLambda > CandidateConfigs.MaxPtLambda) {
+          continue;
+        }
+        if (std::abs(casc.eta()) > CandidateConfigs.etaCasc)
+          continue;
+        if (std::abs(etaLambda) > CandidateConfigs.etaLambdaMax)
+          continue;
+        int bin2DXi = hAcceptanceXi->FindBin(casc.pt(), casc.eta());
+        int bin2DOmega = hAcceptanceOmega->FindBin(casc.pt(), casc.eta());
+        int bin2DLambda = hAcceptanceLambda->FindBin(ptLambda, etaLambda);
+        MeanCos2ThetaLambdaFromXi = hAcceptanceXi->GetBinContent(bin2DXi);
+        MeanCos2ThetaLambdaFromOmega = hAcceptanceOmega->GetBinContent(bin2DOmega);
+        MeanCos2ThetaProtonFromLambda = hAcceptanceLambda->GetBinContent(bin2DLambda);
+      }
+
+      int ChargeIndex = 0;
+      if (casc.sign() > 0)
+        ChargeIndex = 1;
+      double Pzs2Xi = cosThetaStarLambda[0] * std::sin(2 * (casc.phi() - PsiT0C)) / cascadev2::AlphaXi[ChargeIndex] / MeanCos2ThetaLambdaFromXi;
+      double Pzs2Omega = cosThetaStarLambda[1] * std::sin(2 * (casc.phi() - PsiT0C)) / cascadev2::AlphaOmega[ChargeIndex] / MeanCos2ThetaLambdaFromOmega;
+      double Cos2ThetaXi = cosThetaStarLambda[0] * cosThetaStarLambda[0];
+      double Cos2ThetaOmega = cosThetaStarLambda[1] * cosThetaStarLambda[1];
+      double Pzs2LambdaFromCasc = cosThetaStarProton * std::sin(2 * (casc.phi() - PsiT0C)) / cascadev2::AlphaLambda[ChargeIndex] / MeanCos2ThetaProtonFromLambda;
+      double Cos2ThetaLambda = cosThetaStarProton * cosThetaStarProton;
+
+      double CosThetaXiWithAlpha = cosThetaStarLambda[0] / cascadev2::AlphaXi[ChargeIndex];
+      double CosThetaOmegaWithAlpha = cosThetaStarLambda[1] / cascadev2::AlphaOmega[ChargeIndex];
+      double CosThetaProtonWithAlpha = cosThetaStarProton / cascadev2::AlphaLambda[ChargeIndex];
+
+      histos.fill(HIST("hv2CEPvsFT0C"), coll.centFT0C(), v2CEP);
+      histos.fill(HIST("hv2CEPvsv2CSP"), v2CSP, v2CEP);
+      histos.fill(HIST("hCascadePhi"), casc.phi());
+      histos.fill(HIST("hcascminuspsiT0C"), cascminuspsiT0C);
+
+      double values[4]{casc.mXi(), casc.pt(), v2CSP, coll.centFT0C()};
+      if (isSelectedCasc[0]) {
+        cascadev2::hSparseV2C[0]->Fill(values);
+      }
+      if (isSelectedCasc[1]) {
+        values[0] = casc.mOmega();
+        cascadev2::hSparseV2C[0]->Fill(values);
+      }
+
+      float BDTresponse[nParticles]{0.f, 0.f};
+      if (isApplyML) {
+        BDTresponse[0] = bdtScore[0][1];
+        BDTresponse[1] = bdtScore[1][1];
+      }
+
+      if (std::abs(casc.eta()) < CandidateConfigs.etaCasc) {
+        if (fillingConfigs.isFillTHNXi) {
+          if (fillingConfigs.isFillTHN_V2)
+            histos.get<THn>(HIST("hXiV2"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mXi(), BDTresponse[0], v2CEP);
+          if (fillingConfigs.isFillTHN_Pz)
+            histos.get<THn>(HIST("hXiPzs2"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mXi(), BDTresponse[0], Pzs2Xi);
+          if (casc.mLambda() > CandidateConfigs.MinLambdaMass && casc.mLambda() < CandidateConfigs.MaxLambdaMass) {
+            if (fillingConfigs.isFillTHN_PzFromLambda)
+              histos.get<THn>(HIST("hXiPzs2FromLambda"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mXi(), BDTresponse[0], Pzs2LambdaFromCasc);
+          }
+          if (fillingConfigs.isFillTHN_Acc)
+            histos.get<THn>(HIST("hXiCos2Theta"))->Fill(coll.centFT0C(), ChargeIndex, casc.eta(), casc.pt(), casc.mXi(), BDTresponse[0], Cos2ThetaXi);
+          if (fillingConfigs.isFillTHN_AccFromLambdaVsCasc)
+            histos.get<THn>(HIST("hXiCos2ThetaFromLambda"))->Fill(coll.centFT0C(), ChargeIndex, casc.eta(), casc.pt(), casc.mXi(), BDTresponse[0], Cos2ThetaLambda);
+          if (casc.mXi() > CandidateConfigs.MinXiMass && casc.mXi() < CandidateConfigs.MaxXiMass) {
+            if (fillingConfigs.isFillTHN_AccFromLambdaVsLambda)
+              histos.get<THn>(HIST("hXiCos2ThetaFromLambdaL"))->Fill(coll.centFT0C(), ChargeIndex, etaLambda, ptLambda, casc.mLambda(), BDTresponse[0], Cos2ThetaLambda);
+            histos.get<TH1>(HIST("massXi_ProtonAcc"))->Fill(casc.mXi());
+          }
+        }
+        if (fillingConfigs.isFillTHNXi_PzVsPsi) {
+          if (fillingConfigs.isFillTHN_Pz)
+            histos.get<THn>(HIST("hXiPzVsPsi"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mXi(), BDTresponse[0], CosThetaXiWithAlpha, 2 * cascminuspsiT0C);
+          if (fillingConfigs.isFillTHN_PzFromLambda)
+            histos.get<THn>(HIST("hXiPzVsPsiFromLambda"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mXi(), BDTresponse[0], CosThetaProtonWithAlpha, 2 * cascminuspsiT0C);
+          if (casc.mLambda() > CandidateConfigs.MinLambdaMass && casc.mLambda() < CandidateConfigs.MaxLambdaMass) {
+            if (fillingConfigs.isFillTHN_Acc)
+              histos.get<THn>(HIST("hXiCos2ThetaVsPsi"))->Fill(coll.centFT0C(), ChargeIndex, casc.eta(), casc.pt(), casc.mXi(), BDTresponse[0], Cos2ThetaXi, 2 * cascminuspsiT0C);
+          }
+          if (fillingConfigs.isFillTHN_AccFromLambdaVsCasc)
+            histos.get<THn>(HIST("hXiCos2ThetaVsPsiFromLambda"))->Fill(coll.centFT0C(), ChargeIndex, casc.eta(), casc.pt(), casc.mXi(), BDTresponse[0], Cos2ThetaLambda, 2 * cascminuspsiT0C);
+          if (casc.mXi() > CandidateConfigs.MinXiMass && casc.mXi() < CandidateConfigs.MaxXiMass) {
+            if (fillingConfigs.isFillTHN_AccFromLambdaVsLambda)
+              histos.get<THn>(HIST("hXiCos2ThetaVsPsiFromLambdaL"))->Fill(coll.centFT0C(), ChargeIndex, etaLambda, ptLambda, casc.mLambda(), BDTresponse[0], Cos2ThetaLambda, 2 * cascminuspsiT0C);
+            histos.get<TH1>(HIST("massXi_ProtonAcc"))->Fill(casc.mXi());
+          }
+        }
+        if (fillingConfigs.isFillTHNOmega) {
+          if (fillingConfigs.isFillTHN_V2)
+            histos.get<THn>(HIST("hOmegaV2"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mOmega(), BDTresponse[1], v2CEP);
+          if (fillingConfigs.isFillTHN_Pz)
+            histos.get<THn>(HIST("hOmegaPzs2"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mOmega(), BDTresponse[1], Pzs2Omega);
+          if (casc.mLambda() > CandidateConfigs.MinLambdaMass && casc.mLambda() < CandidateConfigs.MaxLambdaMass) {
+            if (fillingConfigs.isFillTHN_PzFromLambda)
+              histos.get<THn>(HIST("hOmegaPzs2FromLambda"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mOmega(), BDTresponse[1], Pzs2LambdaFromCasc);
+          }
+          if (fillingConfigs.isFillTHN_Acc)
+            histos.get<THn>(HIST("hOmegaCos2Theta"))->Fill(coll.centFT0C(), ChargeIndex, casc.eta(), casc.pt(), casc.mOmega(), BDTresponse[1], Cos2ThetaOmega);
+          if (fillingConfigs.isFillTHN_AccFromLambdaVsCasc)
+            histos.get<THn>(HIST("hOmegaCos2ThetaFromLambda"))->Fill(coll.centFT0C(), ChargeIndex, casc.eta(), casc.pt(), casc.mOmega(), BDTresponse[1], Cos2ThetaLambda);
+          if (casc.mOmega() > CandidateConfigs.MinOmegaMass && casc.mOmega() < CandidateConfigs.MaxOmegaMass) {
+            if (fillingConfigs.isFillTHN_AccFromLambdaVsLambda)
+              histos.get<THn>(HIST("hOmegaCos2ThetaFromLambdaL"))->Fill(coll.centFT0C(), ChargeIndex, etaLambda, ptLambda, casc.mLambda(), BDTresponse[1], Cos2ThetaLambda);
+            histos.get<TH1>(HIST("massOmega_ProtonAcc"))->Fill(casc.mOmega());
+          }
+        }
+        if (fillingConfigs.isFillTHNOmega_PzVsPsi) {
+          if (fillingConfigs.isFillTHN_Pz)
+            histos.get<THn>(HIST("hOmegaPzVsPsi"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mOmega(), BDTresponse[1], CosThetaOmegaWithAlpha, 2 * cascminuspsiT0C);
+          if (fillingConfigs.isFillTHN_PzFromLambda)
+            histos.get<THn>(HIST("hOmegaPzVsPsiFromLambda"))->Fill(coll.centFT0C(), ChargeIndex, casc.pt(), casc.mOmega(), BDTresponse[1], CosThetaProtonWithAlpha, 2 * cascminuspsiT0C);
+          if (casc.mLambda() > CandidateConfigs.MinLambdaMass && casc.mLambda() < CandidateConfigs.MaxLambdaMass) {
+            if (fillingConfigs.isFillTHN_Acc)
+              histos.get<THn>(HIST("hOmegaCos2ThetaVsPsi"))->Fill(coll.centFT0C(), ChargeIndex, casc.eta(), casc.pt(), casc.mOmega(), BDTresponse[1], Cos2ThetaOmega, 2 * cascminuspsiT0C);
+          }
+          if (fillingConfigs.isFillTHN_AccFromLambdaVsCasc)
+            histos.get<THn>(HIST("hOmegaCos2ThetaVsPsiFromLambda"))->Fill(coll.centFT0C(), ChargeIndex, casc.eta(), casc.pt(), casc.mOmega(), BDTresponse[1], Cos2ThetaLambda, 2 * cascminuspsiT0C);
+          if (casc.mOmega() > CandidateConfigs.MinOmegaMass && casc.mOmega() < CandidateConfigs.MaxOmegaMass) {
+            if (fillingConfigs.isFillTHN_AccFromLambdaVsLambda)
+              histos.get<THn>(HIST("hOmegaCos2ThetaVsPsiFromLambdaL"))->Fill(coll.centFT0C(), ChargeIndex, etaLambda, ptLambda, casc.mLambda(), BDTresponse[1], Cos2ThetaLambda, 2 * cascminuspsiT0C);
+            histos.get<TH1>(HIST("massOmega_ProtonAcc"))->Fill(casc.mOmega());
+          }
+        }
+      }
+
+      if (isSelectedCasc[0] || isSelectedCasc[1]) {
+        if (fillingConfigs.isFillTree)
+          fillAnalysedTable(coll, hasEventPlane, 0, casc, v2CSP, v2CEP, 0, 0, PsiT0C, BDTresponse[0], BDTresponse[1], 0);
+      }
+    }
+  }
+
+  void processAnalyseDataEPCentralFW(CollEventAndSpecPlaneCentralFW const& coll, CascCandidates const& Cascades, DauTracks const&)
+  {
+
+    if (!AcceptEvent(coll, 1)) {
+      return;
+    }
+
+    // select only events used for the calibration of the event plane
+    if (isGoodEventEP) {
+      if (std::abs(coll.qvecFT0CRe()) > 990 || std::abs(coll.qvecFT0CIm()) > 990 || std::abs(coll.qvecBNegRe()) > 990 || std::abs(coll.qvecBNegIm()) > 990 || std::abs(coll.qvecBPosRe()) > 990 || std::abs(coll.qvecBPosIm()) > 990) {
+        return;
+      }
+    }
+
+    // event has FT0C event plane
     bool hasEventPlane = 0;
-    if (abs(coll.qvecFT0CRe()) > 990 || abs(coll.qvecFT0CIm()) > 990 || abs(coll.qvecBNegRe()) > 990 || abs(coll.qvecBNegIm()) > 990 || abs(coll.qvecBPosRe()) > 990 || abs(coll.qvecBPosIm()) > 990)
+    if (std::abs(coll.qvecFT0CRe()) < 990 && std::abs(coll.qvecFT0CIm()) < 990)
       hasEventPlane = 1;
 
     // event has spectator plane
@@ -1092,9 +1366,9 @@ struct cascadeFlow {
     ROOT::Math::XYZVector spectatorplaneVecZDCA{std::cos(coll.psiZDCA()), std::sin(coll.psiZDCA()), 0}; // eta positive = projectile
     ROOT::Math::XYZVector spectatorplaneVecZDCC{std::cos(coll.psiZDCC()), std::sin(coll.psiZDCC()), 0}; // eta negative = target
 
-    float NormQvT0C = sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C));
-    float NormQvTPCA = sqrt(eventplaneVecTPCA.Dot(eventplaneVecTPCA));
-    float NormQvTPCC = sqrt(eventplaneVecTPCC.Dot(eventplaneVecTPCC));
+    float NormQvT0C = std::sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C));
+    float NormQvTPCA = std::sqrt(eventplaneVecTPCA.Dot(eventplaneVecTPCA));
+    float NormQvTPCC = std::sqrt(eventplaneVecTPCC.Dot(eventplaneVecTPCC));
 
     const float PsiT0C = std::atan2(coll.qvecFT0CIm(), coll.qvecFT0CRe()) * 0.5f;
     histos.fill(HIST("hPsiT0C"), PsiT0C);
@@ -1108,8 +1382,8 @@ struct cascadeFlow {
     resolution.fill(HIST("QVectorsNormTPCAC"), eventplaneVecTPCA.Dot(eventplaneVecTPCC) / (NormQvTPCA * NormQvTPCC), coll.centFT0C());
     resolution.fill(HIST("QVectorsSpecPlane"), spectatorplaneVecZDCC.Dot(spectatorplaneVecZDCA), coll.centFT0C());
 
-    std::vector<float> bdtScore[2];
-    for (auto& casc : Cascades) {
+    std::vector<float> bdtScore[nParticles];
+    for (auto const& casc : Cascades) {
 
       /// Add some minimal cuts for single track variables (min number of TPC clusters)
       auto negExtra = casc.negTrackExtra_as<DauTracks>();
@@ -1121,7 +1395,7 @@ struct cascadeFlow {
       histos.fill(HIST("hCascade"), counter);
 
       // ML selections
-      bool isSelectedCasc[2]{false, false};
+      bool isSelectedCasc[nParticles]{false, false};
 
       std::vector<float> inputFeaturesCasc{casc.cascradius(),
                                            casc.v0radius(),
@@ -1136,7 +1410,7 @@ struct cascadeFlow {
                                            casc.bachBaryonCosPA(),
                                            casc.bachBaryonDCAxyToPV()};
 
-      float massCasc[2]{casc.mXi(), casc.mOmega()};
+      float massCasc[nParticles]{casc.mXi(), casc.mOmega()};
 
       // inv mass loose cut
       if (casc.pt() < CandidateConfigs.MinPt || casc.pt() > CandidateConfigs.MaxPt) {
@@ -1151,7 +1425,7 @@ struct cascadeFlow {
         isSelectedCasc[0] = mlResponseXi.isSelectedMl(inputFeaturesCasc, casc.pt(), bdtScore[0]);
         isSelectedCasc[1] = mlResponseOmega.isSelectedMl(inputFeaturesCasc, casc.pt(), bdtScore[1]);
 
-        for (int iS{0}; iS < 2; ++iS) {
+        for (int iS{0}; iS < nParticles; ++iS) {
           // Fill BDT score histograms before selection
           cascadev2::hSignalScoreBeforeSel[iS]->Fill(bdtScore[0][1]);
           cascadev2::hBkgScoreBeforeSel[iS]->Fill(bdtScore[1][0]);
@@ -1171,12 +1445,12 @@ struct cascadeFlow {
       ROOT::Math::XYZVector cascQvec{std::cos(2 * casc.phi()), std::sin(2 * casc.phi()), 0};
       auto v2CSP = cascQvec.Dot(eventplaneVecT0C);
       auto cascminuspsiT0C = GetPhiInRange(casc.phi() - PsiT0C);
-      auto v2CEP = TMath::Cos(2.0 * cascminuspsiT0C);
+      auto v2CEP = std::cos(2.0 * cascminuspsiT0C);
       ROOT::Math::XYZVector cascUvec{std::cos(casc.phi()), std::sin(casc.phi()), 0};
       auto v1SP_ZDCA = cascUvec.Dot(spectatorplaneVecZDCA);
       auto v1SP_ZDCC = cascUvec.Dot(spectatorplaneVecZDCC);
-      auto v1EP_ZDCA = TMath::Cos(casc.phi() - coll.psiZDCA());
-      auto v1EP_ZDCC = TMath::Cos(casc.phi() - coll.psiZDCC());
+      auto v1EP_ZDCA = std::cos(casc.phi() - coll.psiZDCA());
+      auto v1EP_ZDCC = std::cos(casc.phi() - coll.psiZDCC());
       float v1SP = 0.5 * (v1SP_ZDCA - v1SP_ZDCC);
       float v1EP = 0.5 * (v1EP_ZDCA - v1EP_ZDCC); // same as v1SP
 
@@ -1195,7 +1469,7 @@ struct cascadeFlow {
         cascadev2::hSparseV2C[0]->Fill(values);
       }
 
-      float BDTresponse[2]{0.f, 0.f};
+      float BDTresponse[nParticles]{0.f, 0.f};
       if (isApplyML) {
         BDTresponse[0] = bdtScore[0][1];
         BDTresponse[1] = bdtScore[1][1];
@@ -1224,8 +1498,8 @@ struct cascadeFlow {
     histos.fill(HIST("hEventVertexZ"), coll.posZ());
     histos.fill(HIST("hMultNTracksITSTPCVsCentrality"), coll.centFT0C(), coll.multNTracksITSTPC());
 
-    std::vector<float> bdtScore[2];
-    for (auto& casc : Cascades) {
+    std::vector<float> bdtScore[nParticles];
+    for (auto const& casc : Cascades) {
 
       if (!casc.has_cascMCCore())
         continue;
@@ -1233,8 +1507,8 @@ struct cascadeFlow {
       auto cascMC = casc.cascMCCore_as<soa::Join<aod::CascMCCores, aod::CascMCCollRefs>>();
 
       int pdgCode{cascMC.pdgCode()};
-      if (!(std::abs(pdgCode) == 3312 && std::abs(cascMC.pdgCodeV0()) == 3122 && std::abs(cascMC.pdgCodeBachelor()) == 211)     // Xi
-          && !(std::abs(pdgCode) == 3334 && std::abs(cascMC.pdgCodeV0()) == 3122 && std::abs(cascMC.pdgCodeBachelor()) == 321)) // Omega
+      if (!(std::abs(pdgCode) == PDG_t::kXiMinus && std::abs(cascMC.pdgCodeV0()) == PDG_t::kLambda0 && std::abs(cascMC.pdgCodeBachelor()) == PDG_t::kPiPlus)       // Xi
+          && !(std::abs(pdgCode) == PDG_t::kOmegaMinus && std::abs(cascMC.pdgCodeV0()) == PDG_t::kLambda0 && std::abs(cascMC.pdgCodeBachelor()) == PDG_t::kKPlus)) // Omega
       {
         pdgCode = 0;
       }
@@ -1243,13 +1517,13 @@ struct cascadeFlow {
       float XiY = RecoDecay::y(std::array{casc.px(), casc.py(), casc.pz()}, constants::physics::MassXiMinus);
       float OmegaY = RecoDecay::y(std::array{casc.px(), casc.py(), casc.pz()}, constants::physics::MassOmegaMinus);
       // true reco cascades before applying any selection
-      if (std::abs(pdgCode) == 3312 && std::abs(cascMC.pdgCodeV0()) == 3122 && std::abs(cascMC.pdgCodeBachelor()) == 211) {
+      if (std::abs(pdgCode) == PDG_t::kXiMinus && std::abs(cascMC.pdgCodeV0()) == PDG_t::kLambda0 && std::abs(cascMC.pdgCodeBachelor()) == PDG_t::kPiPlus) {
         histos.fill(HIST("hXiPtvsCent"), coll.centFT0C(), casc.pt());
         if (std::abs(casc.eta()) < 0.8)
           histos.fill(HIST("hXiPtvsCentEta08"), coll.centFT0C(), casc.pt());
         if (std::abs(XiY) < 0.5)
           histos.fill(HIST("hXiPtvsCentY05"), coll.centFT0C(), casc.pt());
-      } else if (std::abs(pdgCode) == 3334 && std::abs(cascMC.pdgCodeV0()) == 3122 && std::abs(cascMC.pdgCodeBachelor()) == 321) {
+      } else if (std::abs(pdgCode) == PDG_t::kOmegaMinus && std::abs(cascMC.pdgCodeV0()) == PDG_t::kLambda0 && std::abs(cascMC.pdgCodeBachelor()) == PDG_t::kKPlus) {
         histos.fill(HIST("hOmegaPtvsCent"), coll.centFT0C(), casc.pt());
         if (std::abs(casc.eta()) < 0.8)
           histos.fill(HIST("hOmegaPtvsCentEta08"), coll.centFT0C(), casc.pt());
@@ -1267,7 +1541,7 @@ struct cascadeFlow {
       histos.fill(HIST("hCascade"), counter);
 
       // ML selections
-      bool isSelectedCasc[2]{false, false};
+      bool isSelectedCasc[nParticles]{false, false};
 
       std::vector<float> inputFeaturesCasc{casc.cascradius(),
                                            casc.v0radius(),
@@ -1282,7 +1556,7 @@ struct cascadeFlow {
                                            casc.bachBaryonCosPA(),
                                            casc.bachBaryonDCAxyToPV()};
 
-      float massCasc[2]{casc.mXi(), casc.mOmega()};
+      float massCasc[nParticles]{casc.mXi(), casc.mOmega()};
 
       if (casc.pt() < CandidateConfigs.MinPt || casc.pt() > CandidateConfigs.MaxPt) {
         continue;
@@ -1296,7 +1570,7 @@ struct cascadeFlow {
         isSelectedCasc[0] = mlResponseXi.isSelectedMl(inputFeaturesCasc, casc.pt(), bdtScore[0]);
         isSelectedCasc[1] = mlResponseOmega.isSelectedMl(inputFeaturesCasc, casc.pt(), bdtScore[1]);
 
-        for (int iS{0}; iS < 2; ++iS) {
+        for (int iS{0}; iS < nParticles; ++iS) {
           // Fill BDT score histograms before selection
           cascadev2::hSignalScoreBeforeSel[iS]->Fill(bdtScore[0][1]);
           cascadev2::hBkgScoreBeforeSel[iS]->Fill(bdtScore[1][0]);
@@ -1315,7 +1589,7 @@ struct cascadeFlow {
 
       histos.fill(HIST("hCascadePhi"), casc.phi());
 
-      float BDTresponse[2]{0.f, 0.f};
+      float BDTresponse[nParticles]{0.f, 0.f};
       const float PsiT0C = 0; // not defined in MC for now
       auto v2CSP = 0;         // not defined in MC for now
       auto v2CEP = 0;         // not defined in MC for now
@@ -1341,7 +1615,7 @@ struct cascadeFlow {
     histosMCGen.fill(HIST("hZvertexGen"), mcCollision.posZ());
     histosMCGen.fill(HIST("hNEventsMC"), 0.5);
     // Generated with accepted z vertex
-    if (TMath::Abs(mcCollision.posZ()) > cutzvertex) {
+    if (std::abs(mcCollision.posZ()) > cutzvertex) {
       return;
     }
     histosMCGen.fill(HIST("hNEventsMC"), 1.5);
@@ -1374,21 +1648,21 @@ struct cascadeFlow {
     histosMCGen.fill(HIST("hNEventsMC"), 5.5);
 
     for (auto const& cascmc : cascMC) {
-      if (TMath::Abs(cascmc.pdgCode()) == 3312)
+      if (std::abs(cascmc.pdgCode()) == PDG_t::kXiMinus)
         histosMCGen.fill(HIST("hNCascGen"), 0.5);
-      else if (TMath::Abs(cascmc.pdgCode()) == 3334)
+      else if (std::abs(cascmc.pdgCode()) == PDG_t::kOmegaMinus)
         histosMCGen.fill(HIST("hNCascGen"), 1.5);
       if (!cascmc.has_straMCCollision())
         continue;
-      if (TMath::Abs(cascmc.pdgCode()) == 3312)
+      if (std::abs(cascmc.pdgCode()) == PDG_t::kXiMinus)
         histosMCGen.fill(HIST("hNCascGen"), 2.5);
-      else if (TMath::Abs(cascmc.pdgCode()) == 3334)
+      else if (std::abs(cascmc.pdgCode()) == PDG_t::kOmegaMinus)
         histosMCGen.fill(HIST("hNCascGen"), 3.5);
       if (!cascmc.isPhysicalPrimary())
         continue;
-      if (TMath::Abs(cascmc.pdgCode()) == 3312)
+      if (std::abs(cascmc.pdgCode()) == PDG_t::kXiMinus)
         histosMCGen.fill(HIST("hNCascGen"), 4.5);
-      else if (TMath::Abs(cascmc.pdgCode()) == 3334)
+      else if (std::abs(cascmc.pdgCode()) == PDG_t::kOmegaMinus)
         histosMCGen.fill(HIST("hNCascGen"), 5.5);
 
       float ptmc = RecoDecay::sqrtSumOfSquares(cascmc.pxMC(), cascmc.pyMC());
@@ -1402,27 +1676,26 @@ struct cascadeFlow {
         theta1 = theta; // 0 < theta1/2 < pi/4 --> 0 < tan (theta1/2) < 1 --> positive eta
       // if pz is negative (i.e. negative rapidity): -pi/2 < theta < 0 --> we need 0 < theta1/2 < pi/2 for the ln to be defined
       else
-        theta1 = TMath::Pi() + theta; // pi/2 < theta1 < pi --> pi/4 < theta1/2 <  pi/2 --> 1 < tan (theta1/2) --> negative eta
+        theta1 = o2::constants::math::PI + theta; // pi/2 < theta1 < pi --> pi/4 < theta1/2 <  pi/2 --> 1 < tan (theta1/2) --> negative eta
 
-      float cascMCeta = -log(std::tan(theta1 / 2));
+      float cascMCeta = -std::log(std::tan(theta1 / 2));
       float cascMCy = 0;
-
-      if (TMath::Abs(cascmc.pdgCode()) == 3312) {
+      if (std::abs(cascmc.pdgCode()) == PDG_t::kXiMinus) {
         cascMCy = RecoDecay::y(std::array{cascmc.pxMC(), cascmc.pyMC(), cascmc.pzMC()}, constants::physics::MassXiMinus);
-        if (TMath::Abs(cascMCeta) < etaCascMCGen) {
+        if (std::abs(cascMCeta) < etaCascMCGen) {
           histosMCGen.fill(HIST("h2DGenXiEta08"), centrality, ptmc);
           histosMCGen.fill(HIST("hNCascGen"), 6.5);
         }
-        if (TMath::Abs(cascMCy) < yCascMCGen)
+        if (std::abs(cascMCy) < yCascMCGen)
           histosMCGen.fill(HIST("h2DGenXiY05"), centrality, ptmc);
         histosMCGen.fill(HIST("hGenXiY"), cascMCy);
-      } else if (TMath::Abs(cascmc.pdgCode() == 3334)) {
+      } else if (std::abs(cascmc.pdgCode()) == PDG_t::kOmegaMinus) {
         cascMCy = RecoDecay::y(std::array{cascmc.pxMC(), cascmc.pyMC(), cascmc.pzMC()}, constants::physics::MassOmegaMinus);
-        if (TMath::Abs(cascMCeta) < etaCascMCGen) {
+        if (std::abs(cascMCeta) < etaCascMCGen) {
           histosMCGen.fill(HIST("h2DGenOmegaEta08"), centrality, ptmc);
           histosMCGen.fill(HIST("hNCascGen"), 7.5);
         }
-        if (TMath::Abs(cascMCy) < yCascMCGen)
+        if (std::abs(cascMCy) < yCascMCGen)
           histosMCGen.fill(HIST("h2DGenOmegaY05"), centrality, ptmc);
         histosMCGen.fill(HIST("hGenOmegaY"), cascMCy);
       }
@@ -1432,6 +1705,7 @@ struct cascadeFlow {
   PROCESS_SWITCH(cascadeFlow, processTrainingBackground, "Process to create the training dataset for the background", true);
   PROCESS_SWITCH(cascadeFlow, processTrainingSignal, "Process to create the training dataset for the signal", false);
   PROCESS_SWITCH(cascadeFlow, processAnalyseData, "Process to apply ML model to the data", false);
+  PROCESS_SWITCH(cascadeFlow, processAnalyseDataEP2CentralFW, "Process to apply ML model to the data - second order event plane calibration from central framework", false);
   PROCESS_SWITCH(cascadeFlow, processAnalyseDataEPCentralFW, "Process to apply ML model to the data - event plane calibration from central framework", false);
   PROCESS_SWITCH(cascadeFlow, processAnalyseMC, "Process to apply ML model to the MC", false);
   PROCESS_SWITCH(cascadeFlow, processMCGen, "Process to store MC generated particles", false);