diff --git a/PWGCF/Flow/TableProducer/zdcQVectors.cxx b/PWGCF/Flow/TableProducer/zdcQVectors.cxx
index 9c1cce63a0b..0aade4713d6 100644
--- a/PWGCF/Flow/TableProducer/zdcQVectors.cxx
+++ b/PWGCF/Flow/TableProducer/zdcQVectors.cxx
@@ -111,7 +111,7 @@ struct ZdcQVectors {
   } EvSel;
 
   struct : ConfigurableGroup {
-    O2_DEFINE_CONFIGURABLE(cfgRecenterForTimestamp, bool, false, "Add 1D recentering for timestamp");
+    O2_DEFINE_CONFIGURABLE(cfgCCDBDir_RecenterForTimestamp, bool, false, "Add 1D recentering for timestamp");
     O2_DEFINE_CONFIGURABLE(cfgCCDBdir_Timestamp, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/Timestamp", "CCDB directory for Timestamp recentering");
   } extraTS;
 
@@ -131,14 +131,12 @@ struct ZdcQVectors {
 
   O2_DEFINE_CONFIGURABLE(cfgVtxZ, float, 10.0f, "Accepted z-vertex range")
   O2_DEFINE_CONFIGURABLE(cfgMagField, float, 99999, "Configurable magnetic field; default CCDB will be queried")
-  O2_DEFINE_CONFIGURABLE(cfgEnergyCal, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/Energy", "ccdb path for energy calibration histos")
-  O2_DEFINE_CONFIGURABLE(cfgMeanv, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/vmean", "ccdb path for mean v histos")
   O2_DEFINE_CONFIGURABLE(cfgMinEntriesSparseBin, int, 1000, "Minimal number of entries allowed in 4D recentering histogram to use for recentering.")
-  O2_DEFINE_CONFIGURABLE(cfgRec, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5", "ccdb path for recentering histos");
   O2_DEFINE_CONFIGURABLE(cfgFillHistRegistry, bool, true, "Fill common registry with histograms");
   O2_DEFINE_CONFIGURABLE(cfgFillCutAnalysis, bool, true, "Fill cut analysis with histograms");
   O2_DEFINE_CONFIGURABLE(cfgFillNothing, bool, false, "Disable ALL Histograms -> ONLY use to reduce memory");
   O2_DEFINE_CONFIGURABLE(cfgNoGain, bool, true, "Do not apply gain correction to ZDC energy calibration");
+  O2_DEFINE_CONFIGURABLE(cfgNIterationsAfterShift, int, 1, "Number of iterations to perform after shift");
 
   O2_DEFINE_CONFIGURABLE(cfgTrackSelsDCAxy, float, 0.2, "Cut on DCA in the transverse direction (cm)");
   O2_DEFINE_CONFIGURABLE(cfgTrackSelsDCAz, float, 0.2, "Cut on DCA in the longitudinal direction (cm)");
@@ -146,11 +144,19 @@ struct ZdcQVectors {
   O2_DEFINE_CONFIGURABLE(cfgTrackSelsPtmax, float, 10, "maximum pt (GeV/c)");
   O2_DEFINE_CONFIGURABLE(cfgTrackSelsEta, float, 0.8, "eta cut");
 
+  O2_DEFINE_CONFIGURABLE(cfgCCDBDir_EnergyCal, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/Energy", "ccdb path for energy calibration histos")
+  O2_DEFINE_CONFIGURABLE(cfgCCDBDir_Meanv, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/vmean", "ccdb path for mean v histos")
+  O2_DEFINE_CONFIGURABLE(cfgCCDBDir_Rec, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5", "ccdb path for recentering histos");
   O2_DEFINE_CONFIGURABLE(cfgCCDBdir_Shift, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/Shift", "CCDB directory for Shift ZDC");
+  O2_DEFINE_CONFIGURABLE(cfgCCDBdir_ShiftRec, std::string, "Users/c/ckoster/ZDC/LHC23_PbPb_pass5/noGain/AfterShift", "CCDB directory for recentering after Shift ZDC");
+
   Configurable<std::vector<int>> cfgSelVec{"cfgSelVec", std::vector<int>{1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1}, "Put 1 for every event selection from SelectionCriteria that is used in flowSP"};
   Configurable<std::vector<double>> cfgEvSelsMultPv{"cfgEvSelsMultPv", std::vector<double>{2223.49, -75.1444, 0.963572, -0.00570399, 1.34877e-05, 3790.99, -137.064, 2.13044, -0.017122, 5.82834e-05}, "Multiplicity cuts (PV) first 5 parameters cutLOW last 5 cutHIGH (Default is +-2sigma pass5) "};
   Configurable<std::vector<double>> cfgEvSelsMult{"cfgEvSelsMult", std::vector<double>{1301.56, -41.4615, 0.478224, -0.00239449, 4.46966e-06, 2967.6, -102.927, 1.47488, -0.0106534, 3.28622e-05}, "Multiplicity cuts (Global) first 5 parameters cutLOW last 5 cutHIGH (Default is +-2sigma pass5) "};
 
+  // Track selection DCA cut
+  std::unique_ptr<TF1> fDCACut = std::make_unique<TF1>("fDCACut", "0.0105 + 0.035 / TMath::Power(x,1.1)", 0, 100);
+
   // define my.....
   // Filter collisionFilter = nabs(aod::collision::posZ) <;
 
@@ -185,6 +191,7 @@ struct ZdcQVectors {
     kMeanv,
     kRec,
     kTimestamp,
+    kRecShift,
     nCalibModes
   };
 
@@ -204,10 +211,8 @@ struct ZdcQVectors {
 
   // keep track of calibration histos for each given step and iteration
   struct Calib {
-    std::vector<TList*> calibList = std::vector<TList*>(4, nullptr); // [0] Enerfy cal, [1] vmean, [2] recentering, [3] timestamp
-    std::vector<bool> calibfilesLoaded = std::vector<bool>(4, false);
-    int atStep = 0;
-    int atIteration = 0;
+    std::vector<TList*> calibList = std::vector<TList*>(5, nullptr); // [0] Enerfy cal, [1] vmean, [2] recentering, [3] timestamp, [4] recenter after shift Correction
+    std::vector<bool> calibfilesLoaded = std::vector<bool>(5, false);
 
     TProfile3D* shiftprofileC = nullptr;
     TProfile3D* shiftprofileA = nullptr;
@@ -220,11 +225,13 @@ struct ZdcQVectors {
     uint64_t timestamp = 0;
     std::vector<float> v = {0, 0, 0};
     bool isSelected = 0;
+    int atIteration = 0;
   } cal;
 
   enum FillType {
     kBefore,
-    kAfter
+    kAfter,
+    kAfterShift
   };
 
   void init(InitContext const&)
@@ -270,6 +277,9 @@ struct ZdcQVectors {
         registry.add<TH2>(Form("QA/before/hSPplaneA"), "hSPplaneA", kTH2D, {axisPsiA, axisCent});
         registry.add<TH2>(Form("QA/before/hSPplaneC"), "hSPplaneC", kTH2D, {axisPsiC, axisCent});
         registry.add<TH2>(Form("QA/before/hSPplaneFull"), "hSPplaneFull", kTH2D, {{100, -PI, PI}, axisCent});
+        if (!cfgCCDBdir_ShiftRec.value.empty()) {
+          registry.addClone("QA/before/", "QA/afterShift/");
+        }
         for (const auto& side : sides) {
           registry.add<TH2>(Form("recentering/before/hZN%s_Qx_vs_Qy", side), Form("hZN%s_Qx_vs_Qy", side), kTH2F, {axisQ, axisQ});
         }
@@ -294,7 +304,7 @@ struct ZdcQVectors {
             registry.add(Form("recentering/before/hQ%s%s_vs_vy", coord, side), Form("hQ%s%s_vs_vy", coord, side), {HistType::kTProfile, {axisVy}});
             registry.add(Form("recentering/before/hQ%s%s_vs_vz", coord, side), Form("hQ%s%s_vs_vz", coord, side), {HistType::kTProfile, {axisVz}});
             registry.add(Form("recentering/before/hQ%s%s_vs_timestamp", coord, side), Form("hQ%s%s_vs_timestamp", coord, side), {HistType::kTProfile, {axisTimestamp}});
-            registry.add(Form("recentering/Q%s%s_vs_iteration", coord, side), Form("hQ%s%s_vs_iteration", coord, side), {HistType::kTH2D, {{35, 0, 35}, axisQ}});
+            registry.add(Form("recentering/Q%s%s_vs_iteration", coord, side), Form("hQ%s%s_vs_iteration", coord, side), {HistType::kTH2D, {{100, 0, 100}, axisQ}});
           } // end of capCOORDS
         } // end of sides
 
@@ -377,6 +387,10 @@ struct ZdcQVectors {
         registry.add<TProfile2D>("CutAnalysis/hvertex_vy", "hvertex_vy", kTProfile2D, {{1, 0., 1.}, {nEventSelections + 5, 0, nEventSelections + 5}});
         registry.add<TProfile2D>("CutAnalysis/hvertex_vz", "hvertex_vz", kTProfile2D, {{1, 0., 1.}, {nEventSelections + 5, 0, nEventSelections + 5}});
       }
+
+      if (!cfgCCDBdir_ShiftRec.value.empty()) {
+        registry.addClone("recentering/before/", "recentering/afterShift/");
+      }
       registry.addClone("recentering/before/", "recentering/after/");
       registry.addClone("QA/before/", "QA/after/");
     }
@@ -614,7 +628,7 @@ struct ZdcQVectors {
     if (cfgFillNothing) {
       return;
     }
-    static constexpr std::array<std::string_view, 2> Time = {"before", "after"}; // todo move to struct like in flowSP
+    static constexpr std::array<std::string_view, 3> Time = {"before", "after", "afterShift"};
 
     registry.fill(HIST("recentering/") + HIST(Time[ft]) + HIST("/hZNA_Qx_vs_Qy"), qxa, qya);
     registry.fill(HIST("recentering/") + HIST(Time[ft]) + HIST("/hZNC_Qx_vs_Qy"), qxc, qyc);
@@ -675,11 +689,11 @@ struct ZdcQVectors {
     registry.fill(HIST("recentering/") + HIST(Time[ft]) + HIST("/ZNC_Qy_vs_Centrality"), centrality, qyc);
 
     // add psi!!
-    double psiA = 1.0 * std::atan2(qxc, qxa);
+    double psiA = 1.0 * std::atan2(qya, qxa);
     registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/hSPplaneA"), psiA, centrality, 1);
-    double psiC = 1.0 * std::atan2(qyc, qya);
+    double psiC = 1.0 * std::atan2(qyc, qxc);
     registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/hSPplaneC"), psiC, centrality, 1);
-    double psiFull = 1.0 * std::atan2(qxc + qyc, qxa + qya);
+    double psiFull = 1.0 * std::atan2(qya + qyc, qxa + qxc);
     registry.fill(HIST("QA/") + HIST(Time[ft]) + HIST("/hSPplaneFull"), psiFull, centrality, 1);
   }
 
@@ -698,7 +712,6 @@ struct ZdcQVectors {
       cal.calibfilesLoaded[cm] = true;
       LOGF(info, "Loaded calibration histos from %s", ccdb_dir.c_str());
       if (cm == kRec) {
-        cal.atStep = 5;
         cal.atIteration = 5;
       }
     }
@@ -722,7 +735,7 @@ struct ZdcQVectors {
       if (!hist) {
         LOGF(fatal, "No calibration histo for iteration %i and step %i -> %s", iteration, step, objName);
       }
-    } else if (cm == kRec) {
+    } else if (cm == kRec || cm == kRecShift) {
       auto list = dynamic_cast<TList*>(cal.calibList[cm]->FindObject(Form("it%i_step%i", iteration, step)));
       if (!list) {
         LOGF(fatal, "No calibration list for iteration %i and step %i", iteration, step);
@@ -731,8 +744,6 @@ struct ZdcQVectors {
       if (!hist) {
         LOGF(fatal, "No calibration histo for iteration %i and step %i -> %s", iteration, step, objName);
       }
-      cal.atStep = step;
-      cal.atIteration = iteration;
     }
     if (!hist) {
       LOGF(fatal, "%s not available.. Abort..", objName);
@@ -768,7 +779,7 @@ struct ZdcQVectors {
           bin = h->GetXaxis()->FindBin(TString::Format("%i", cal.runnumber));
         }
         if (name.Contains("timestamp")) {
-          bin = h->GetXaxis()->FindBin(cal.timestamp);
+          bin = h->GetXaxis()->FindBin(cal.rsTimestamp);
         }
         calibConstant = h->GetBinContent(bin);
       }
@@ -817,7 +828,6 @@ struct ZdcQVectors {
     std::vector<float> cents;
 
     auto cent = collision.centFT0C();
-
     cents.push_back(collision.centFT0C());
 
     if (cfgFT0Cvariant1) {
@@ -853,6 +863,30 @@ struct ZdcQVectors {
     cal.runnumber = runnumber;
     cal.centrality = cent;
 
+    // load new calibrations for new runs only
+    if (runnumber != cal.lastRunNumber) {
+      cal.calibfilesLoaded[kEnergyCal] = false;
+      cal.calibList[kEnergyCal] = nullptr;
+
+      cal.calibfilesLoaded[kMeanv] = false;
+      cal.calibList[kMeanv] = nullptr;
+
+      cal.calibfilesLoaded[kRec] = false;
+      cal.calibList[kRec] = nullptr;
+
+      cal.calibfilesLoaded[kTimestamp] = false;
+      cal.calibList[kTimestamp] = nullptr;
+
+      cal.calibfilesLoaded[kRecShift] = false;
+      cal.calibList[kRecShift] = nullptr;
+
+      cal.isShiftProfileFound = false;
+      cal.shiftprofileC = nullptr;
+      cal.shiftprofileA = nullptr;
+
+      cal.atIteration = 0;
+    }
+
     if (cfgFillHistRegistry && !cfgFillNothing) {
       registry.fill(HIST("QA/centrality_before"), cent);
     }
@@ -891,16 +925,20 @@ struct ZdcQVectors {
     bool isZNChit = true;
 
     for (int i = 0; i < nTowers; ++i) {
-      if (i < nTowersPerSide && eZN[i] <= 0)
+      if (i < nTowersPerSide && eZN[i] <= 0) {
         isZNAhit = false;
-      if (i >= nTowersPerSide && eZN[i] <= 0)
+      }
+      if (i >= nTowersPerSide && eZN[i] <= 0) {
         isZNChit = false;
+      }
     }
 
-    if (zdcCol.energyCommonZNA() <= 0)
+    if (zdcCol.energyCommonZNA() <= 0) {
       isZNAhit = false;
-    if (zdcCol.energyCommonZNC() <= 0)
+    }
+    if (zdcCol.energyCommonZNC() <= 0) {
       isZNChit = false;
+    }
 
     // if ZNA or ZNC not hit correctly.. do not use event in q-vector calculation
     if (!isZNAhit || !isZNChit) {
@@ -927,32 +965,13 @@ struct ZdcQVectors {
     }
     registry.fill(HIST("hEventCount"), evSel_CentCuts);
 
-    // load new calibrations for new runs only
-    if (runnumber != cal.lastRunNumber) {
-      cal.calibfilesLoaded[kEnergyCal] = false;
-      cal.calibList[kEnergyCal] = nullptr;
-
-      cal.calibfilesLoaded[kMeanv] = false;
-      cal.calibList[kMeanv] = nullptr;
-
-      cal.calibfilesLoaded[kRec] = false;
-      cal.calibList[kRec] = nullptr;
-
-      cal.calibfilesLoaded[kTimestamp] = false;
-      cal.calibList[kTimestamp] = nullptr;
-
-      cal.isShiftProfileFound = false;
-      cal.shiftprofileC = nullptr;
-      cal.shiftprofileA = nullptr;
-    }
-
     // load the calibration histos for iteration 0 step 0 (Energy Calibration)
     if (!cfgNoGain) {
-      loadCalibrations<kEnergyCal>(cfgEnergyCal.value, timestamp);
+      loadCalibrations<kEnergyCal>(cfgCCDBDir_EnergyCal.value, timestamp);
     }
 
     // load the calibrations for the mean v
-    loadCalibrations<kMeanv>(cfgMeanv.value, timestamp);
+    loadCalibrations<kMeanv>(cfgCCDBDir_Meanv.value, timestamp);
 
     if (!cfgFillNothing && isEventSelected) {
       registry.get<TProfile>(HIST("vmean/hvertex_vx"))->Fill(Form("%d", runnumber), v[0]);
@@ -1074,17 +1093,25 @@ struct ZdcQVectors {
       return;
     }
 
-    loadCalibrations<kRec>(cfgRec.value, timestamp);
+    // Load in TList with all recentering histos
+    loadCalibrations<kRec>(cfgCCDBDir_Rec.value, timestamp);
 
-    if (extraTS.cfgRecenterForTimestamp) {
+    // If CCDB dir given load recentering for extra step Timestamp
+    if (extraTS.cfgCCDBDir_RecenterForTimestamp) {
       loadCalibrations<kTimestamp>(extraTS.cfgCCDBdir_Timestamp.value, timestamp);
     }
 
+    // If CCDB dir given load recentering after shift correction.
+    if (!cfgCCDBdir_ShiftRec.value.empty()) {
+      loadCalibrations<kRecShift>(cfgCCDBdir_ShiftRec.value, timestamp);
+    }
+
     std::array<double, 4> qRec(q);
 
     if (cal.atIteration == 0) {
-      if (cal.isSelected && cfgFillHistRegistry && isEventSelected)
+      if (cal.isSelected && cfgFillHistRegistry && isEventSelected) {
         fillCommonRegistry<kBefore>(q[0], q[1], q[2], q[3], cal.v, cent, rsTimestamp);
+      }
 
       spTableZDC(runnumber, cents, cal.v, foundBC.timestamp(), q[0], q[1], q[2], q[3], cal.isSelected, eventSelectionFlags);
       cal.lastRunNumber = runnumber;
@@ -1136,7 +1163,7 @@ struct ZdcQVectors {
           pb++;
         }
 
-        if (extraTS.cfgRecenterForTimestamp) {
+        if (extraTS.cfgCCDBDir_RecenterForTimestamp) {
           corrQxA.push_back(getCorrection<TProfile, kTimestamp>(namesTS[0].Data(), it, 6));
           corrQyA.push_back(getCorrection<TProfile, kTimestamp>(namesTS[1].Data(), it, 6));
           corrQxC.push_back(getCorrection<TProfile, kTimestamp>(namesTS[2].Data(), it, 6));
@@ -1257,20 +1284,87 @@ struct ZdcQVectors {
       double qXcShift = std::hypot(qRec[2], qRec[3]) * std::cos(psiZDCCshift);
       double qYcShift = std::hypot(qRec[2], qRec[3]) * std::sin(psiZDCCshift);
 
+      if (!cal.calibfilesLoaded[kRecShift]) {
+        if (cal.isSelected && cfgFillHistRegistry && !cfgFillNothing && isEventSelected) {
+          fillCommonRegistry<kAfter>(qXaShift, qYaShift, qXcShift, qYcShift, cal.v, cent, rsTimestamp);
+          registry.fill(HIST("QA/centrality_after"), cent);
+          registry.get<TProfile>(HIST("QA/after/ZNA_Qx"))->Fill(Form("%d", runnumber), qXaShift);
+          registry.get<TProfile>(HIST("QA/after/ZNA_Qy"))->Fill(Form("%d", runnumber), qYaShift);
+          registry.get<TProfile>(HIST("QA/after/ZNC_Qx"))->Fill(Form("%d", runnumber), qXcShift);
+          registry.get<TProfile>(HIST("QA/after/ZNC_Qy"))->Fill(Form("%d", runnumber), qYcShift);
+        }
+
+        spTableZDC(runnumber, cents, cal.v, foundBC.timestamp(), qXaShift, qYaShift, qXcShift, qYcShift, cal.isSelected, eventSelectionFlags);
+        qRec = {0, 0, 0, 0};
+
+        cal.lastRunNumber = runnumber;
+        return;
+      }
+      // vector of 4
+      corrQxA.clear();
+      corrQyA.clear();
+      corrQxC.clear();
+      corrQyC.clear();
+
       if (cal.isSelected && cfgFillHistRegistry && !cfgFillNothing && isEventSelected) {
         fillCommonRegistry<kAfter>(qXaShift, qYaShift, qXcShift, qYcShift, cal.v, cent, rsTimestamp);
-        registry.fill(HIST("QA/centrality_after"), cent);
-        registry.get<TProfile>(HIST("QA/after/ZNA_Qx"))->Fill(Form("%d", runnumber), qXaShift);
-        registry.get<TProfile>(HIST("QA/after/ZNA_Qy"))->Fill(Form("%d", runnumber), qYaShift);
-        registry.get<TProfile>(HIST("QA/after/ZNC_Qx"))->Fill(Form("%d", runnumber), qXcShift);
-        registry.get<TProfile>(HIST("QA/after/ZNC_Qy"))->Fill(Form("%d", runnumber), qYcShift);
       }
 
-      spTableZDC(runnumber, cents, cal.v, foundBC.timestamp(), qXaShift, qYaShift, qXcShift, qYcShift, cal.isSelected, eventSelectionFlags);
-      qRec = {0, 0, 0, 0};
+      for (int it = 1; it <= cfgNIterationsAfterShift; it++) {
+        corrQxA.push_back(getCorrection<THnSparse, kRecShift>(names[0][0].Data(), it, 1));
+        corrQyA.push_back(getCorrection<THnSparse, kRecShift>(names[0][1].Data(), it, 1));
+        corrQxC.push_back(getCorrection<THnSparse, kRecShift>(names[0][2].Data(), it, 1));
+        corrQyC.push_back(getCorrection<THnSparse, kRecShift>(names[0][3].Data(), it, 1));
 
-      cal.lastRunNumber = runnumber;
-      return;
+        if (cfgFillHistRegistry && !cfgFillNothing && isEventSelected) {
+          registry.get<TH2>(HIST("recentering/QXA_vs_iteration"))->Fill(pb + 1, qXaShift - std::accumulate(corrQxA.begin(), corrQxA.end(), 0.0));
+          registry.get<TH2>(HIST("recentering/QYA_vs_iteration"))->Fill(pb + 1, qYaShift - std::accumulate(corrQyA.begin(), corrQyA.end(), 0.0));
+          registry.get<TH2>(HIST("recentering/QXC_vs_iteration"))->Fill(pb + 1, qXcShift - std::accumulate(corrQxC.begin(), corrQxC.end(), 0.0));
+          registry.get<TH2>(HIST("recentering/QYC_vs_iteration"))->Fill(pb + 1, qYcShift - std::accumulate(corrQyC.begin(), corrQyC.end(), 0.0));
+        }
+        pb++;
+
+        for (int step = 2; step <= nSteps; step++) {
+          corrQxA.push_back(getCorrection<TProfile, kRecShift>(names[step - 1][0].Data(), it, step));
+          corrQyA.push_back(getCorrection<TProfile, kRecShift>(names[step - 1][1].Data(), it, step));
+          corrQxC.push_back(getCorrection<TProfile, kRecShift>(names[step - 1][2].Data(), it, step));
+          corrQyC.push_back(getCorrection<TProfile, kRecShift>(names[step - 1][3].Data(), it, step));
+
+          if (cfgFillHistRegistry && !cfgFillNothing && isEventSelected) {
+            registry.get<TH2>(HIST("recentering/QXA_vs_iteration"))->Fill(pb + 1, q[0] - std::accumulate(corrQxA.begin(), corrQxA.end(), 0.0));
+            registry.get<TH2>(HIST("recentering/QYA_vs_iteration"))->Fill(pb + 1, q[1] - std::accumulate(corrQyA.begin(), corrQyA.end(), 0.0));
+            registry.get<TH2>(HIST("recentering/QXC_vs_iteration"))->Fill(pb + 1, q[2] - std::accumulate(corrQxC.begin(), corrQxC.end(), 0.0));
+            registry.get<TH2>(HIST("recentering/QYC_vs_iteration"))->Fill(pb + 1, q[3] - std::accumulate(corrQyC.begin(), corrQyC.end(), 0.0));
+          }
+
+          pb++;
+        }
+
+        double totalCorrectionQxAshift = std::accumulate(corrQxA.begin(), corrQxA.end(), 0.0);
+        double totalCorrectionQyAshift = std::accumulate(corrQyA.begin(), corrQyA.end(), 0.0);
+        double totalCorrectionQxCshift = std::accumulate(corrQxC.begin(), corrQxC.end(), 0.0);
+        double totalCorrectionQyCshift = std::accumulate(corrQyC.begin(), corrQyC.end(), 0.0);
+
+        qXaShift -= totalCorrectionQxAshift;
+        qYaShift -= totalCorrectionQyAshift;
+        qXcShift -= totalCorrectionQxCshift;
+        qYcShift -= totalCorrectionQyCshift;
+
+        if (cal.isSelected && cfgFillHistRegistry && !cfgFillNothing && isEventSelected) {
+          fillCommonRegistry<kAfterShift>(qXaShift, qYaShift, qXcShift, qYcShift, cal.v, cent, rsTimestamp);
+          registry.fill(HIST("QA/centrality_after"), cent);
+          registry.get<TProfile>(HIST("QA/after/ZNA_Qx"))->Fill(Form("%d", runnumber), qXaShift);
+          registry.get<TProfile>(HIST("QA/after/ZNA_Qy"))->Fill(Form("%d", runnumber), qYaShift);
+          registry.get<TProfile>(HIST("QA/after/ZNC_Qx"))->Fill(Form("%d", runnumber), qXcShift);
+          registry.get<TProfile>(HIST("QA/after/ZNC_Qy"))->Fill(Form("%d", runnumber), qYcShift);
+        }
+
+        spTableZDC(runnumber, cents, cal.v, foundBC.timestamp(), qXaShift, qYaShift, qXcShift, qYcShift, cal.isSelected, eventSelectionFlags);
+        qRec = {0, 0, 0, 0};
+
+        cal.lastRunNumber = runnumber;
+        return;
+      }
 
     LOGF(warning, "We return without saving table... -> THis is a problem");
     cal.lastRunNumber = runnumber;
diff --git a/PWGCF/Flow/Tasks/flowSP.cxx b/PWGCF/Flow/Tasks/flowSP.cxx
index 9e04b351488..26b96ec9c80 100644
--- a/PWGCF/Flow/Tasks/flowSP.cxx
+++ b/PWGCF/Flow/Tasks/flowSP.cxx
@@ -53,6 +53,7 @@
 #include <TProfile2D.h>
 #include <TProfile3D.h>
 
+#include <array>
 #include <chrono>
 #include <cmath>
 #include <cstdint>
@@ -71,7 +72,7 @@ using namespace o2::framework::expressions;
 using namespace o2::aod::rctsel;
 // using namespace o2::analysis;
 
-#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};
+#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP}; // NOLINT(bugprone-macro-parentheses)
 
 struct FlowSP {
 
@@ -130,6 +131,7 @@ struct FlowSP {
     O2_DEFINE_CONFIGURABLE(cCentMin, float, 0, "Minimum cenrality for selected events");
     O2_DEFINE_CONFIGURABLE(cCentMax, float, 90, "Maximum cenrality for selected events");
     O2_DEFINE_CONFIGURABLE(cFilterLeptons, bool, true, "Filter out leptons from MCGenerated by requiring |pdgCode| > 100");
+    O2_DEFINE_CONFIGURABLE(cDoGeneratedInReco, bool, true, "Do MC Generated inside the MCReco process function on mcCollisions");
     // NUA and NUE weights
     O2_DEFINE_CONFIGURABLE(cFillWeights, bool, true, "Fill NUA weights");
     O2_DEFINE_CONFIGURABLE(cFillWeightsPOS, bool, true, "Fill NUA weights only for positive charges");
@@ -142,9 +144,7 @@ struct FlowSP {
     // Additional track Selections
     O2_DEFINE_CONFIGURABLE(cTrackSelsUseAdditionalTrackCut, bool, false, "Bool to enable Additional Track Cut");
     O2_DEFINE_CONFIGURABLE(cTrackSelsDoDCApt, bool, true, "Apply Pt dependent DCAz cut");
-    O2_DEFINE_CONFIGURABLE(cTrackSelsDCApt1, float, 0.1, "DcaZ < const + (a * b) / pt^1.1 -> this sets a");
-    O2_DEFINE_CONFIGURABLE(cTrackSelsDCApt2, float, 0.035, "DcaZ < const + (a * b) / pt^1.1 -> this sets b");
-    O2_DEFINE_CONFIGURABLE(cTrackSelsDCAptConsMin, float, 0.1, "DcaZ < const + (a * b) / pt^1.1 -> this sets const");
+    O2_DEFINE_CONFIGURABLE(cTrackSelsDCAfunc, std::string, "0.0105 + 0.035 / TMath::Power(x,1.1)", "Function for pT-dependent DCA-cut (Default y = 0.0105 + 0.035 / TMath::Power(x,1.1))");
     O2_DEFINE_CONFIGURABLE(cTrackSelsPIDNsigma, float, 2.0, "nSigma cut for PID");
     O2_DEFINE_CONFIGURABLE(cTrackSelDoTrackQAvsCent, bool, true, "Do track selection QA plots as function of centrality");
     // harmonics for v coefficients
@@ -176,8 +176,8 @@ struct FlowSP {
   Filter trackFilter = nabs(aod::track::eta) < cfg.cTrackSelsEta && aod::track::pt > cfg.cTrackSelsPtmin&& aod::track::pt < cfg.cTrackSelsPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true) || cfg.cIsMCReco) && nabs(aod::track::dcaXY) < cfg.cTrackSelsDCAxy&& nabs(aod::track::dcaZ) < cfg.cTrackSelsDCAz;
   Filter trackFilterMC = nabs(aod::mcparticle::eta) < cfg.cTrackSelsEta && aod::mcparticle::pt > cfg.cTrackSelsPtmin&& aod::mcparticle::pt < cfg.cTrackSelsPtmax;
   using GeneralCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::CentNGlobals>;
-  using UnfilteredTracksPID = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTOFbeta, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr>;
-  using UnfilteredTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA>;
+  using UnfilteredTracksPID = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTOFbeta, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr>;
+  using UnfilteredTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA, aod::TrackSelectionExtension>;
 
   using UsedTracks = soa::Filtered<UnfilteredTracks>;
   using UsedTracksPID = soa::Filtered<UnfilteredTracksPID>;
@@ -197,16 +197,16 @@ struct FlowSP {
   Preslice<TCs> trackPerCollision = aod::track::collisionId;
 
   //  Connect to ccdb
-  Service<ccdb::BasicCCDBManager> ccdb;
-  Service<o2::framework::O2DatabasePDG> pdg;
+  Service<ccdb::BasicCCDBManager> ccdb{};
+  Service<o2::framework::O2DatabasePDG> pdg{};
 
   // struct to hold the correction histos/
   struct Config {
-    std::vector<TH1D*> mEfficiency = {};
-    std::vector<TH2D*> mEfficiency2D = {};
-    std::vector<TH3D*> mEfficiency3D = {};
-    std::vector<GFWWeights*> mAcceptance = {};
-    std::vector<TH3D*> mAcceptance2D = {};
+    std::vector<TH1D*> mEfficiency;
+    std::vector<TH2D*> mEfficiency2D;
+    std::vector<TH3D*> mEfficiency3D;
+    std::vector<GFWWeights*> mAcceptance;
+    std::vector<TH3D*> mAcceptance2D;
     bool correctionsLoaded = false;
     int lastRunNumber = 0;
 
@@ -290,6 +290,9 @@ struct FlowSP {
   std::unique_ptr<TF1> fMultCutHigh = nullptr;
   std::unique_ptr<TF1> fMultMultPVCut = nullptr;
 
+  // Track selection DCA cut
+  std::unique_ptr<TF1> fDCACut = std::make_unique<TF1>("fDCACut", cfg.cTrackSelsDCAfunc.value.c_str(), 0, 100);
+
   enum SelectionCriteria {
     evSel_FilteredEvent,
     evSel_sel8,
@@ -312,12 +315,12 @@ struct FlowSP {
     trackSel_ZeroCharge,
     trackSel_Eta,
     trackSel_Pt,
+    trackSel_NCls,
+    trackSel_TPCBoundary,
     trackSel_DCAxy,
     trackSel_DCAz,
-    trackSel_GlobalTracks,
-    trackSel_NCls,
     trackSel_FshCls,
-    trackSel_TPCBoundary,
+    trackSel_GlobalTracks,
     trackSel_ParticleWeights,
     nTrackSelections
   };
@@ -349,9 +352,9 @@ struct FlowSP {
     nParticleTypes
   };
 
-  static constexpr std::string_view Charge[] = {"incl/", "pos/", "neg/"};
-  static constexpr std::string_view Species[] = {"", "pion/", "kaon/", "proton/"};
-  static constexpr std::string_view Time[] = {"before/", "after/"};
+  static constexpr std::array<std::string_view, 3> Charge = {"incl/", "pos/", "neg/"};
+  static constexpr std::array<std::string_view, 4> Species = {"", "pion/", "kaon/", "proton/"};
+  static constexpr std::array<std::string_view, 2> Time = {"before/", "after/"};
 
   void init(InitContext const&)
   {
@@ -499,8 +502,9 @@ struct FlowSP {
           histos.add<TH3>("incl/QA/after/hPhi_Eta_Pt_corrected", "", kTH3D, {axisPhi, axisEta, axisPt});
         }
 
-        if (cfg.cFillQABefore)
+        if (cfg.cFillQABefore) {
           histos.addClone("incl/QA/after/", "incl/QA/before/");
+        }
       }
 
       if (cfg.cFillPIDQA && doprocessDataPID) {
@@ -548,9 +552,22 @@ struct FlowSP {
       }
 
       if (doprocessMCReco) {
-        registry.add("trackMCReco/after/incl/hIsPhysicalPrimary", "", {HistType::kTH3D, {{2, 0, 2}, axisCentrality, axisPt}});
-        registry.get<TH3>(HIST("trackMCReco/after/incl/hIsPhysicalPrimary"))->GetXaxis()->SetBinLabel(1, "Secondary");
-        registry.get<TH3>(HIST("trackMCReco/after/incl/hIsPhysicalPrimary"))->GetXaxis()->SetBinLabel(2, "Primary");
+        registry.add("trackMCReco/after/incl/hIsPhysicalPrimary", "", {HistType::kTH3D, {axisCentrality, axisEta, axisPt}});
+        registry.add("trackMCReco/after/incl/hIsNotPhysicalPrimary", "", {HistType::kTH3D, {axisCentrality, axisEta, axisPt}});
+        registry.add("trackMCReco/incl/hPtMCPtTrack", "", {HistType::kTH2D, {axisPt, axisPt}});
+        registry.add("trackMCReco/incl/hEtaMCEtaTrack", "", {HistType::kTH2D, {axisEta, axisEta}});
+        registry.add("trackMCReco/incl/hPtPerTrackSelection", "", {HistType::kTH2D, {axisPt, {nTrackSelections, 0, nTrackSelections}}});
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_Eta + 1, "Eta");
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_Pt + 1, "Pt");
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_DCAxy + 1, "DCAxy");
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_DCAz + 1, "DCAz");
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_GlobalTracks + 1, "GlobalTracks");
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_NCls + 1, "nClusters TPC");
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_FshCls + 1, "Frac. sh. Cls TPC");
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_TPCBoundary + 1, "TPC Boundary");
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_ZeroCharge + 1, "Only charged");
+        registry.get<TH2>(HIST("trackMCReco/incl/hPtPerTrackSelection"))->GetYaxis()->SetBinLabel(trackSel_ParticleWeights + 1, "Apply weights");
+
         registry.add("trackMCReco/hTrackSize_unFiltered", "", {HistType::kTH2D, {{100, 0, 4000}, axisCentrality}});
         registry.add("trackMCReco/hTrackSize_Filtered", "", {HistType::kTH2D, {{100, 0, 4000}, axisCentrality}});
         registry.add("trackMCReco/after/incl/hPt_hadron", "", {HistType::kTH3D, {axisPt, axisEta, axisCentrality}});
@@ -560,6 +577,8 @@ struct FlowSP {
         // Clone into particles and before/after
         registry.addClone("trackMCReco/after/incl/", "trackMCReco/after/pos/");
         registry.addClone("trackMCReco/after/incl/", "trackMCReco/after/neg/");
+        registry.addClone("trackMCReco/incl/", "trackMCReco/neg/");
+        registry.addClone("trackMCReco/incl/", "trackMCReco/pos/");
         registry.addClone("trackMCReco/after/", "trackMCReco/before/");
       }
 
@@ -681,14 +700,15 @@ struct FlowSP {
           registry.addClone("incl/", "neg/");
         }
       }
-
-    } else if (doprocessMCGen) {
+    }
+    if (doprocessMCGen || doprocessMCReco) {
       registry.add("trackMCGen/nCollReconstructedPerMcCollision", "", {HistType::kTH1D, {{10, -5, 5}}});
       registry.add("trackMCGen/after/incl/hPt_hadron", "", {HistType::kTH3D, {axisPt, axisEta, axisCentrality}});
       registry.add("trackMCGen/after/incl/hPt_proton", "", {HistType::kTH3D, {axisPt, axisEta, axisCentrality}});
       registry.add("trackMCGen/after/incl/hPt_pion", "", {HistType::kTH3D, {axisPt, axisEta, axisCentrality}});
       registry.add("trackMCGen/after/incl/hPt_kaon", "", {HistType::kTH3D, {axisPt, axisEta, axisCentrality}});
       registry.add("trackMCGen/after/incl/phi_eta_vtxZ_gen", "", {HistType::kTH3D, {axisPhi, axisEta, axisVz}});
+
       registry.addClone("trackMCGen/after/incl/", "trackMCGen/after/pos/");
       registry.addClone("trackMCGen/after/incl/", "trackMCGen/after/neg/");
       registry.addClone("trackMCGen/after/", "trackMCGen/before/");
@@ -768,8 +788,9 @@ struct FlowSP {
     int etaind = hNUA->GetYaxis()->FindBin(eta);
     int vzind = hNUA->GetZaxis()->FindBin(vtxz);
     float weight = hNUA->GetBinContent(xind, etaind, vzind);
-    if (weight != 0)
+    if (weight != 0) {
       return 1. / weight;
+    }
     return 1;
   }
 
@@ -809,13 +830,11 @@ struct FlowSP {
       }
     }
 
-    if (nIdentified == 0) {
-      return kUnidentified; // No PID match found
-    } else if (nIdentified == 1) {
+    if (nIdentified == 1) {
       return valPID;
-    } else {
-      return kUnidentified; // Multiple PID matches found
     }
+
+    return kUnidentified; // Multiple PID matches found
   }
 
   int getMagneticField(uint64_t timestamp)
@@ -850,28 +869,30 @@ struct FlowSP {
   void loadCorrections(uint64_t timestamp)
   {
     // corrections saved on CCDB as TList {incl, pos, neg} of GFWWeights (acc) TH1D (eff) objects!
-    if (conf.correctionsLoaded)
+    if (conf.correctionsLoaded) {
       return;
+    }
 
     int nWeights = 3;
 
     if (cfg.cUseNUA1D) {
-      if (cfg.cCCDB_NUA.value.empty() == false) {
-        TList* listCorrections = ccdb->getForTimeStamp<TList>(cfg.cCCDB_NUA, timestamp);
-        conf.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights")));
-        conf.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights_positive")));
-        conf.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights_negative")));
+      if (!cfg.cCCDB_NUA.value.empty()) {
+        auto* listCorrections = ccdb->getForTimeStamp<TList>(cfg.cCCDB_NUA, timestamp);
+        conf.mAcceptance.push_back(dynamic_cast<GFWWeights*>(listCorrections->FindObject("weights")));
+        conf.mAcceptance.push_back(dynamic_cast<GFWWeights*>(listCorrections->FindObject("weights_positive")));
+        conf.mAcceptance.push_back(dynamic_cast<GFWWeights*>(listCorrections->FindObject("weights_negative")));
         int sizeAcc = conf.mAcceptance.size();
-        if (sizeAcc < nWeights)
+        if (sizeAcc < nWeights) {
           LOGF(fatal, "Could not load acceptance weights from %s", cfg.cCCDB_NUA.value.c_str());
-        else
+        } else {
           LOGF(info, "Loaded acceptance weights from %s", cfg.cCCDB_NUA.value.c_str());
+        }
       } else {
         LOGF(info, "cfg.cCCDB_NUA empty! No corrections loaded");
       }
     } else if (cfg.cUseNUA2D) {
-      if (cfg.cCCDB_NUA.value.empty() == false) {
-        TH3D* hNUA2D = ccdb->getForTimeStamp<TH3D>(cfg.cCCDB_NUA, timestamp);
+      if (!cfg.cCCDB_NUA.value.empty()) {
+        auto* hNUA2D = ccdb->getForTimeStamp<TH3D>(cfg.cCCDB_NUA, timestamp);
         if (!hNUA2D) {
           LOGF(fatal, "Could not load acceptance weights from %s", cfg.cCCDB_NUA.value.c_str());
         } else {
@@ -883,39 +904,41 @@ struct FlowSP {
       }
     }
     // Get Efficiency correction
-    if (cfg.cCCDB_NUE.value.empty() == false) {
-      TList* listCorrections = ccdb->getForTimeStamp<TList>(cfg.cCCDB_NUE, timestamp);
-      conf.mEfficiency.push_back(reinterpret_cast<TH1D*>(listCorrections->FindObject("Efficiency")));
-      conf.mEfficiency.push_back(reinterpret_cast<TH1D*>(listCorrections->FindObject("Efficiency_pos")));
-      conf.mEfficiency.push_back(reinterpret_cast<TH1D*>(listCorrections->FindObject("Efficiency_neg")));
+    if (!cfg.cCCDB_NUE.value.empty()) {
+      auto* listCorrections = ccdb->getForTimeStamp<TList>(cfg.cCCDB_NUE, timestamp);
+      conf.mEfficiency.push_back(dynamic_cast<TH1D*>(listCorrections->FindObject("Efficiency")));
+      conf.mEfficiency.push_back(dynamic_cast<TH1D*>(listCorrections->FindObject("Efficiency_pos")));
+      conf.mEfficiency.push_back(dynamic_cast<TH1D*>(listCorrections->FindObject("Efficiency_neg")));
       int sizeEff = conf.mEfficiency.size();
-      if (sizeEff < nWeights)
+      if (sizeEff < nWeights) {
         LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfg.cCCDB_NUE.value.c_str());
-      else
+      } else {
         LOGF(info, "Loaded efficiency histogram from %s", cfg.cCCDB_NUE.value.c_str());
+      }
     } else {
       LOGF(info, "cfg.cCCDB_NUE empty! No corrections loaded");
     }
     // Get Efficiency correction
-    if (cfg.cCCDB_NUE2D.value.empty() == false) {
-      TList* listCorrections = ccdb->getForTimeStamp<TList>(cfg.cCCDB_NUE2D, timestamp);
-      conf.mEfficiency2D.push_back(reinterpret_cast<TH2D*>(listCorrections->FindObject("Efficiency")));
-      conf.mEfficiency2D.push_back(reinterpret_cast<TH2D*>(listCorrections->FindObject("Efficiency_pos")));
-      conf.mEfficiency2D.push_back(reinterpret_cast<TH2D*>(listCorrections->FindObject("Efficiency_neg")));
+    if (!cfg.cCCDB_NUE2D.value.empty()) {
+      auto* listCorrections = ccdb->getForTimeStamp<TList>(cfg.cCCDB_NUE2D, timestamp);
+      conf.mEfficiency2D.push_back(dynamic_cast<TH2D*>(listCorrections->FindObject("Efficiency")));
+      conf.mEfficiency2D.push_back(dynamic_cast<TH2D*>(listCorrections->FindObject("Efficiency_pos")));
+      conf.mEfficiency2D.push_back(dynamic_cast<TH2D*>(listCorrections->FindObject("Efficiency_neg")));
       int sizeEff = conf.mEfficiency2D.size();
-      if (sizeEff < nWeights)
+      if (sizeEff < nWeights) {
         LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfg.cCCDB_NUE.value.c_str());
-      else
+      } else {
         LOGF(info, "Loaded efficiency histogram from %s", cfg.cCCDB_NUE.value.c_str());
+      }
     } else {
       LOGF(info, "cfg.cCCDB_NUE2 empty! No corrections loaded");
     }
 
-    if (cfg.cCCDB_NUE3D.value.empty() == false) {
-      TList* listCorrections = ccdb->getForTimeStamp<TList>(cfg.cCCDB_NUE3D, timestamp);
-      conf.mEfficiency3D.push_back(reinterpret_cast<TH3D*>(listCorrections->FindObject("Efficiency")));
-      conf.mEfficiency3D.push_back(reinterpret_cast<TH3D*>(listCorrections->FindObject("Efficiency_pos")));
-      conf.mEfficiency3D.push_back(reinterpret_cast<TH3D*>(listCorrections->FindObject("Efficiency_neg")));
+    if (!cfg.cCCDB_NUE3D.value.empty()) {
+      auto* listCorrections = ccdb->getForTimeStamp<TList>(cfg.cCCDB_NUE3D, timestamp);
+      conf.mEfficiency3D.push_back(dynamic_cast<TH3D*>(listCorrections->FindObject("Efficiency")));
+      conf.mEfficiency3D.push_back(dynamic_cast<TH3D*>(listCorrections->FindObject("Efficiency_pos")));
+      conf.mEfficiency3D.push_back(dynamic_cast<TH3D*>(listCorrections->FindObject("Efficiency_neg")));
       int sizeEff = conf.mEfficiency3D.size();
       if (sizeEff < nWeights)
         LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfg.cCCDB_NUE.value.c_str());
@@ -977,15 +1000,28 @@ struct FlowSP {
     return true;
   }
 
+  template <FillType ft, ChargeType ct, typename McParticleObject>
+  inline void fillPrimaryHistos(const McParticleObject& mcparticle)
+  {
+
+    if (mcparticle.isPhysicalPrimary()) {
+      registry.fill(HIST("trackMCReco/") + HIST(Time[ft]) + HIST(Charge[ct]) + HIST("hIsPhysicalPrimary"), spm.centrality, mcparticle.eta(), mcparticle.pt());
+    } else {
+      registry.fill(HIST("trackMCReco/") + HIST(Time[ft]) + HIST(Charge[ct]) + HIST("hIsNotPhysicalPrimary"), spm.centrality, mcparticle.eta(), mcparticle.pt());
+    }
+  }
+
   template <typename TCollision>
   bool eventSelected(const TCollision& collision, const int& multTrk)
   {
-    if (!collision.sel8())
+    if (!collision.sel8()) {
       return 0;
+    }
     histos.fill(HIST("hEventCount"), evSel_sel8);
 
-    if (cfg.cEvtUseRCTFlagChecker && !rctChecker(collision))
+    if (cfg.cEvtUseRCTFlagChecker && !rctChecker(collision)) {
       return 0;
+    }
     histos.fill(HIST("hEventCount"), evSel_RCTFlagsZDC);
 
     // Occupancy
@@ -994,8 +1030,8 @@ struct FlowSP {
       if (occupancy > cfg.cEvSelsMaxOccupancy || occupancy < cfg.cEvSelsMinOccupancy) {
         return 0;
       }
-      histos.fill(HIST("hEventCount"), evSel_occupancy);
     }
+    histos.fill(HIST("hEventCount"), evSel_occupancy);
 
     if (cfg.cEvSelsNoSameBunchPileupCut) {
       if (!collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) {
@@ -1003,37 +1039,37 @@ struct FlowSP {
         // https://indico.cern.ch/event/1396220/#1-event-selection-with-its-rof
         return 0;
       }
-      histos.fill(HIST("hEventCount"), evSel_kNoSameBunchPileup);
     }
+    histos.fill(HIST("hEventCount"), evSel_kNoSameBunchPileup);
     if (cfg.cEvSelsIsGoodZvtxFT0vsPV) {
       if (!collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) {
         // removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference
         // use this cut at low multiplicities with caution
         return 0;
       }
-      histos.fill(HIST("hEventCount"), evSel_kIsGoodZvtxFT0vsPV);
     }
+    histos.fill(HIST("hEventCount"), evSel_kIsGoodZvtxFT0vsPV);
     if (cfg.cEvSelsNoCollInTimeRangeStandard) {
       if (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) {
         //  Rejection of the collisions which have other events nearby
         return 0;
       }
-      histos.fill(HIST("hEventCount"), evSel_kNoCollInTimeRangeStandard);
     }
+    histos.fill(HIST("hEventCount"), evSel_kNoCollInTimeRangeStandard);
     if (cfg.cEvSelsNoCollInTimeRangeNarrow) {
       if (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeNarrow)) {
         // Rejection of the collisions which have other events nearby
         return 0;
       }
-      histos.fill(HIST("hEventCount"), evSel_kNoCollInTimeRangeNarrow);
     }
+    histos.fill(HIST("hEventCount"), evSel_kNoCollInTimeRangeNarrow);
     if (cfg.cEvSelsIsVertexITSTPC) {
       if (!collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) {
         // selects collisions with at least one ITS-TPC track, and thus rejects vertices built from ITS-only tracks
         return 0;
       }
-      histos.fill(HIST("hEventCount"), evSel_kIsVertexITSTPC);
     }
+    histos.fill(HIST("hEventCount"), evSel_kIsVertexITSTPC);
 
     if (cfg.cEvSelsIsGoodITSLayersAll) {
       if (!collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) {
@@ -1041,14 +1077,14 @@ struct FlowSP {
         // https://indico.cern.ch/event/1493023/ (09-01-2025)
         return 0;
       }
-      histos.fill(HIST("hEventCount"), evSel_kIsGoodITSLayersAll);
     }
+    histos.fill(HIST("hEventCount"), evSel_kIsGoodITSLayersAll);
     if (cfg.cEvSelsIsGoodITSLayer0123) {
       if (!collision.selection_bit(o2::aod::evsel::kIsGoodITSLayer0123)) {
         return 0;
       }
-      histos.fill(HIST("hEventCount"), evSel_kIsGoodITSLayer0123);
     }
+    histos.fill(HIST("hEventCount"), evSel_kIsGoodITSLayer0123);
 
     if (cfg.cEvSelsUseAdditionalEventCut) {
       float vtxz = -999;
@@ -1057,25 +1093,30 @@ struct FlowSP {
         float zRes = std::sqrt(collision.covZZ());
         float minzRes = 0.25;
         int maxNumContrib = 20;
-        if (zRes > minzRes && collision.numContrib() < maxNumContrib)
+        if (zRes > minzRes && collision.numContrib() < maxNumContrib) {
           vtxz = -999;
+        }
       }
 
       auto multNTracksPV = collision.multNTracksPV();
 
-      if (vtxz > cfg.cEvSelsVtxZ || vtxz < -cfg.cEvSelsVtxZ)
+      if (vtxz > cfg.cEvSelsVtxZ || vtxz < -cfg.cEvSelsVtxZ) {
         return 0;
-      if (multNTracksPV < fMultPVCutLow->Eval(collision.centFT0C()))
+      }
+      if (multNTracksPV < fMultPVCutLow->Eval(collision.centFT0C())) {
         return 0;
-      if (multNTracksPV > fMultPVCutHigh->Eval(collision.centFT0C()))
+      }
+      if (multNTracksPV > fMultPVCutHigh->Eval(collision.centFT0C())) {
         return 0;
-      if (multTrk < fMultCutLow->Eval(collision.centFT0C()))
+      }
+      if (multTrk < fMultCutLow->Eval(collision.centFT0C())) {
         return 0;
-      if (multTrk > fMultCutHigh->Eval(collision.centFT0C()))
+      }
+      if (multTrk > fMultCutHigh->Eval(collision.centFT0C())) {
         return 0;
-
-      histos.fill(HIST("hEventCount"), evSel_MultCuts);
+      }
     }
+    histos.fill(HIST("hEventCount"), evSel_MultCuts);
 
     return 1;
   }
@@ -1083,64 +1124,131 @@ struct FlowSP {
   template <typename TrackObject>
   bool trackSelected(const TrackObject& track, const int& field)
   {
-    if (std::fabs(track.eta()) > cfg.cTrackSelsEta)
+    auto fillSpectraStudyMCReco = [&](TrackSelections sel) {
+      if constexpr (o2::framework::has_type_v<aod::mctracklabel::McParticleId, typename TrackObject::all_columns>) {
+        if (!track.has_mcParticle()) {
+          return;
+        }
+        const auto mcParticle = track.template mcParticle_as<MCs>();
+        if (doprocessMCReco && mcParticle.isPhysicalPrimary()) {
+          registry.fill(HIST("trackMCReco/incl/hPtPerTrackSelection"), mcParticle.pt(), sel);
+          if (track.sign() > 0) {
+            registry.fill(HIST("trackMCReco/pos/hPtPerTrackSelection"), mcParticle.pt(), sel);
+          } else if (track.sign() < 0) {
+            registry.fill(HIST("trackMCReco/neg/hPtPerTrackSelection"), mcParticle.pt(), sel);
+          }
+        }
+      }
+    };
+
+    if (std::fabs(track.eta()) > cfg.cTrackSelsEta) {
       return false;
+    }
     histos.fill(HIST("hTrackCount"), trackSel_Eta);
+    fillSpectraStudyMCReco(trackSel_Eta);
 
-    if (track.pt() < cfg.cTrackSelsPtmin || track.pt() > cfg.cTrackSelsPtmax)
+    if (track.pt() < cfg.cTrackSelsPtmin || track.pt() > cfg.cTrackSelsPtmax) {
       return false;
-
+    }
     histos.fill(HIST("hTrackCount"), trackSel_Pt);
+    fillSpectraStudyMCReco(trackSel_Pt);
 
-    if (track.dcaXY() > cfg.cTrackSelsDCAxy)
-      return false;
-
-    histos.fill(HIST("hTrackCount"), trackSel_DCAxy);
-
-    if (track.dcaZ() > cfg.cTrackSelsDCAz)
-      return false;
-
-    if (cfg.cTrackSelsDoDCApt && std::fabs(track.dcaZ()) > (cfg.cTrackSelsDCAptConsMin + (cfg.cTrackSelsDCApt1 * cfg.cTrackSelsDCApt2) / (std::pow(track.pt(), 1.1))))
-      return false;
-
-    histos.fill(HIST("hTrackCount"), trackSel_DCAz);
-
-    if (track.tpcNClsCrossedRows() < cfg.cTrackSelsNcls)
+    if (track.tpcNClsCrossedRows() < cfg.cTrackSelsNcls) {
       return false;
+    }
     histos.fill(HIST("hTrackCount"), trackSel_NCls);
-
-    if (track.tpcFractionSharedCls() > cfg.cTrackSelsFshcls)
-      return false;
-    histos.fill(HIST("hTrackCount"), trackSel_FshCls);
+    fillSpectraStudyMCReco(trackSel_NCls);
 
     double phimodn = track.phi();
-    if (field < 0) // for negative polarity field
+    if (field < 0) {
       phimodn = o2::constants::math::TwoPI - phimodn;
-    if (track.sign() < 0) // for negative charge
+    }
+    if (track.sign() < 0) {
       phimodn = o2::constants::math::TwoPI - phimodn;
-    if (phimodn < 0)
+    }
+    if (phimodn < 0) {
       LOGF(warning, "phi < 0: %g", phimodn);
+    }
 
     phimodn += o2::constants::math::PI / 18.0; // to center gap in the middle
     phimodn = fmod(phimodn, o2::constants::math::PI / 9.0);
-    if (cfg.cFillTrackQA && cfg.cFillQABefore)
+    if (cfg.cFillTrackQA && cfg.cFillQABefore) {
       histos.fill(HIST("incl/QA/before/pt_phi"), track.pt(), phimodn);
+    }
 
     if (cfg.cTrackSelsUseAdditionalTrackCut) {
-      if (phimodn < fPhiCutHigh->Eval(track.pt()) && phimodn > fPhiCutLow->Eval(track.pt()))
+      if (phimodn < fPhiCutHigh->Eval(track.pt()) && phimodn > fPhiCutLow->Eval(track.pt())) {
         return false; // reject track
+      }
     }
-    if (cfg.cFillTrackQA)
+    if (cfg.cFillTrackQA) {
       histos.fill(HIST("incl/QA/after/pt_phi"), track.pt(), phimodn);
+    }
     histos.fill(HIST("hTrackCount"), trackSel_TPCBoundary);
+    fillSpectraStudyMCReco(trackSel_TPCBoundary);
+
+    // Only fill primary/secondary histos for MC data.
+    if constexpr (o2::framework::has_type_v<aod::mctracklabel::McParticleId, typename TrackObject::all_columns>) {
+      auto mcParticle = track.template mcParticle_as<MCs>();
+
+      fillPrimaryHistos<kBefore, kInclusive>(mcParticle);
+      if (spm.charge == kPositive) {
+        fillPrimaryHistos<kBefore, kPositive>(mcParticle);
+      } else {
+        fillPrimaryHistos<kBefore, kNegative>(mcParticle);
+      }
+
+      if (std::fabs(track.dcaXY()) > cfg.cTrackSelsDCAxy || std::fabs(track.dcaXY()) > fDCACut->Eval(track.pt())) {
+        return false;
+      }
+      histos.fill(HIST("hTrackCount"), trackSel_DCAxy);
+      fillSpectraStudyMCReco(trackSel_DCAxy);
+
+      if (std::fabs(track.dcaZ()) > cfg.cTrackSelsDCAz || (cfg.cTrackSelsDoDCApt && std::fabs(track.dcaZ()) > fDCACut->Eval(track.pt()))) {
+        return false;
+      }
+      histos.fill(HIST("hTrackCount"), trackSel_DCAz);
+      fillSpectraStudyMCReco(trackSel_DCAz);
+
+      fillPrimaryHistos<kAfter, kInclusive>(mcParticle);
+      if (spm.charge == kPositive) {
+        fillPrimaryHistos<kAfter, kPositive>(mcParticle);
+      } else {
+        fillPrimaryHistos<kAfter, kNegative>(mcParticle);
+      }
+    } else { // Only apply DCA-cuts for data
+      if (std::fabs(track.dcaXY()) > cfg.cTrackSelsDCAxy || std::fabs(track.dcaXY()) > fDCACut->Eval(track.pt())) {
+        return false;
+      }
+      histos.fill(HIST("hTrackCount"), trackSel_DCAxy);
+
+      if (std::fabs(track.dcaZ()) > cfg.cTrackSelsDCAz || (cfg.cTrackSelsDoDCApt && std::fabs(track.dcaZ()) > fDCACut->Eval(track.pt()))) {
+        return false;
+      }
+      histos.fill(HIST("hTrackCount"), trackSel_DCAz);
+    }
+
+    if (track.tpcFractionSharedCls() > cfg.cTrackSelsFshcls) {
+      return false;
+    }
+    histos.fill(HIST("hTrackCount"), trackSel_FshCls);
+    fillSpectraStudyMCReco(trackSel_FshCls);
+
+    if (!track.isGlobalTrack()) {
+      return false;
+    }
+    histos.fill(HIST("hTrackCount"), trackSel_GlobalTracks);
+    fillSpectraStudyMCReco(trackSel_GlobalTracks);
+
     return true;
   }
 
   template <FillType ft, typename CollisionObject, typename TracksObject>
   inline void fillEventQA(const CollisionObject& collision, const TracksObject& tracks)
   {
-    if (!cfg.cFillEventQA)
+    if (!cfg.cFillEventQA) {
       return;
+    }
 
     histos.fill(HIST("QA/") + HIST(Time[ft]) + HIST("hCentFT0C"), collision.centFT0C(), spm.centWeight);
     histos.fill(HIST("QA/") + HIST(Time[ft]) + HIST("hCentNGlobal"), collision.centNGlobal(), spm.centWeight);
@@ -1253,8 +1361,9 @@ struct FlowSP {
   template <FillType ft, ChargeType ct, ParticleType par, typename TrackObject>
   inline void fillTrackQA(const TrackObject& track)
   {
-    if (!cfg.cFillTrackQA)
+    if (!cfg.cFillTrackQA) {
       return;
+    }
 
     double weight = spm.wacc[ct][par] * spm.weff[ct][par] * spm.centWeight;
 
@@ -1282,8 +1391,9 @@ struct FlowSP {
   template <FillType ft, ChargeType ct, typename TrackObject>
   inline void fillPIDQA(const TrackObject& track)
   {
-    if (!cfg.cFillTrackQA)
+    if (!cfg.cFillTrackQA) {
       return;
+    }
     if constexpr (framework::has_type_v<aod::pidtof::TOFNSigmaPi, typename TrackObject::all_columns>) {
       histos.fill(HIST(Charge[ct]) + HIST("pion/") + HIST("QA/") + HIST(Time[ft]) + HIST("hNsigmaTOF_pt"), track.pt(), track.tofNSigmaPi());
       histos.fill(HIST(Charge[ct]) + HIST("pion/") + HIST("QA/") + HIST(Time[ft]) + HIST("hNsigmaTPC_pt"), track.pt(), track.tpcNSigmaPi());
@@ -1332,17 +1442,6 @@ struct FlowSP {
     }
   }
 
-  template <FillType ft, ChargeType ct, typename McParticleObject>
-  inline void fillPrimaryHistos(const McParticleObject& mcparticle)
-  {
-
-    if (!mcparticle.isPhysicalPrimary()) {
-      registry.fill(HIST("trackMCReco/") + HIST(Time[ft]) + HIST(Charge[ct]) + HIST("hIsPhysicalPrimary"), 0, spm.centrality, mcparticle.pt());
-    } else {
-      registry.fill(HIST("trackMCReco/") + HIST(Time[ft]) + HIST(Charge[ct]) + HIST("hIsPhysicalPrimary"), 1, spm.centrality, mcparticle.pt());
-    }
-  }
-
   template <FillType ft, ParticleType par, typename TrackObject>
   void fillAllQA(const TrackObject& track)
   {
@@ -1381,27 +1480,34 @@ struct FlowSP {
       LOGF(info, "Size of mAcceptance: %i (should be 0)", (int)conf.mAcceptance.size());
     }
 
-    if (cfg.cFillQABefore)
+    if (cfg.cFillQABefore) {
       fillEventQA<kBefore>(collision, tracks);
+    }
 
     loadCorrections(bc.timestamp());
 
     spm.centrality = collision.centFT0C();
 
-    if (cfg.cCentFT0Cvariant1)
+    if (cfg.cCentFT0Cvariant1) {
       spm.centrality = collision.centFT0CVariant1();
-    if (cfg.cCentFT0M)
+    }
+    if (cfg.cCentFT0M) {
       spm.centrality = collision.centFT0M();
-    if (cfg.cCentFV0A)
+    }
+    if (cfg.cCentFV0A) {
       spm.centrality = collision.centFV0A();
-    if (cfg.cCentNGlobal)
+    }
+    if (cfg.cCentNGlobal) {
       spm.centrality = collision.centNGlobal();
+    }
 
-    if (!eventSelected(collision, tracks.size()))
+    if (!eventSelected(collision, tracks.size())) {
       return;
+    }
 
-    if (!collision.isSelected()) // selected by ZDCQVectors task (checks signal in ZDC) --> only possible in data not MC
+    if (!collision.isSelected()) { // selected by ZDCQVectors task (checks signal in ZDC) --> only possible in data not MC
       return;
+    }
     histos.fill(HIST("hEventCount"), evSel_isSelectedZDC);
 
     // Always fill centrality histogram after event selections!
@@ -1450,8 +1556,9 @@ struct FlowSP {
       histos.fill(HIST("QA/hFullEvPlaneRes"), spm.centrality, -1 * std::cos(spm.psiA - spm.psiC));
     }
 
-    if (spm.centrality > cfg.cCentMax || spm.centrality < cfg.cCentMin)
+    if (spm.centrality > cfg.cCentMax || spm.centrality < cfg.cCentMin) {
       return;
+    }
 
     histos.fill(HIST("hEventCount"), evSel_CentCuts);
 
@@ -1459,12 +1566,12 @@ struct FlowSP {
     // Only load once!
     // If not loaded set to 1
 
-    if (cfg.cCCDBdir_QQ.value.empty() == false) {
+    if (!cfg.cCCDBdir_QQ.value.empty()) {
       if (!conf.clQQ) {
-        TList* hcorrList = ccdb->getForTimeStamp<TList>(cfg.cCCDBdir_QQ.value, bc.timestamp());
-        conf.hcorrQQ = reinterpret_cast<TProfile*>(hcorrList->FindObject("qAqCXY"));
-        conf.hcorrQQx = reinterpret_cast<TProfile*>(hcorrList->FindObject("qAqCX"));
-        conf.hcorrQQy = reinterpret_cast<TProfile*>(hcorrList->FindObject("qAqCY"));
+        auto* hcorrList = ccdb->getForTimeStamp<TList>(cfg.cCCDBdir_QQ.value, bc.timestamp());
+        conf.hcorrQQ = dynamic_cast<TProfile*>(hcorrList->FindObject("qAqCXY"));
+        conf.hcorrQQx = dynamic_cast<TProfile*>(hcorrList->FindObject("qAqCX"));
+        conf.hcorrQQy = dynamic_cast<TProfile*>(hcorrList->FindObject("qAqCY"));
         conf.clQQ = true;
       }
       spm.corrQQ = conf.hcorrQQ->GetBinContent(conf.hcorrQQ->FindBin(spm.centrality));
@@ -1473,19 +1580,20 @@ struct FlowSP {
     }
 
     double evPlaneRes = 1.;
-    if (cfg.cCCDBdir_SP.value.empty() == false) {
+    if (!cfg.cCCDBdir_SP.value.empty()) {
       if (!conf.clEvPlaneRes) {
         conf.hEvPlaneRes = ccdb->getForTimeStamp<TProfile>(cfg.cCCDBdir_SP.value, bc.timestamp());
         conf.clEvPlaneRes = true;
       }
       evPlaneRes = conf.hEvPlaneRes->GetBinContent(conf.hEvPlaneRes->FindBin(spm.centrality));
-      if (evPlaneRes < 0)
+      if (evPlaneRes < 0) {
         LOGF(fatal, "<Cos(PsiA-PsiC)> > 0 for centrality %.2f! Cannot determine resolution.. Change centrality ranges!!!", spm.centrality);
+      }
       evPlaneRes = std::sqrt(evPlaneRes);
     }
 
     spm.centWeight = 1.;
-    if (cfg.cCCDBdir_centrality.value.empty() == false) {
+    if (!cfg.cCCDBdir_centrality.value.empty()) {
       if (!conf.clCentrality) {
         conf.hCentrality = ccdb->getForTimeStamp<TH1D>(cfg.cCCDBdir_centrality.value, bc.timestamp());
         conf.clCentrality = true;
@@ -1523,8 +1631,9 @@ struct FlowSP {
 
     for (const auto& track : tracks) {
 
-      if (track.sign() == 0)
+      if (track.sign() == 0) {
         continue;
+      }
 
       histos.fill(HIST("hTrackCount"), trackSel_ZeroCharge);
 
@@ -1534,11 +1643,12 @@ struct FlowSP {
         fillAllQA<kBefore, kUnidentified>(track);
       }
 
-      if (!trackSelected(track, field))
+      if (!trackSelected(track, field)) {
         continue;
+      }
 
       spm.meanPtWeight = 1.0;
-      if (cfg.cCCDBdir_meanPt.value.empty() == false) {
+      if (!cfg.cCCDBdir_meanPt.value.empty()) {
         if (!conf.clMeanPt) {
           conf.hMeanPt = ccdb->getForTimeStamp<TProfile2D>(cfg.cCCDBdir_meanPt.value, bc.timestamp());
           conf.clMeanPt = true;
@@ -1573,10 +1683,12 @@ struct FlowSP {
       }
 
       // Set weff and wacc for inclusive, negative and positive hadrons
-      if (!setCurrentParticleWeights(kInclusive, kUnidentified, phi, track.eta(), track.pt(), vtxz, spm.centrality))
+      if (!setCurrentParticleWeights(kInclusive, kUnidentified, phi, track.eta(), track.pt(), vtxz, spm.centrality)) {
         continue;
-      if (!setCurrentParticleWeights(spm.charge, kUnidentified, phi, track.eta(), track.pt(), vtxz, spm.centrality))
+      }
+      if (!setCurrentParticleWeights(spm.charge, kUnidentified, phi, track.eta(), track.pt(), vtxz, spm.centrality)) {
         continue;
+      }
 
       histos.fill(HIST("hTrackCount"), trackSel_ParticleWeights);
 
@@ -1713,20 +1825,26 @@ struct FlowSP {
 
     spm.centrality = collision.centFT0C();
 
-    if (cfg.cCentFT0Cvariant1)
+    if (cfg.cCentFT0Cvariant1) {
       spm.centrality = collision.centFT0CVariant1();
-    if (cfg.cCentFT0M)
+    }
+    if (cfg.cCentFT0M) {
       spm.centrality = collision.centFT0M();
-    if (cfg.cCentFV0A)
+    }
+    if (cfg.cCentFV0A) {
       spm.centrality = collision.centFV0A();
-    if (cfg.cCentNGlobal)
+    }
+    if (cfg.cCentNGlobal) {
       spm.centrality = collision.centNGlobal();
+    }
 
-    if (!eventSelected(collision, tracks.size()))
+    if (!eventSelected(collision, tracks.size())) {
       return;
+    }
 
-    if (!collision.isSelected()) // selected by ZDCQVectors task (checks signal in ZDC) --> only possible in data not MC
+    if (!collision.isSelected()) { // selected by ZDCQVectors task (checks signal in ZDC) --> only possible in data not MC
       return;
+    }
     histos.fill(HIST("hEventCount"), evSel_isSelectedZDC);
 
     // Always fill centrality histogram after event selections!
@@ -1745,19 +1863,20 @@ struct FlowSP {
     // https://twiki.cern.ch/twiki/pub/ALICE/DirectedFlowAnalysisNote/vn_ZDC_ALICE_INT_NOTE_version02.pdf
     spm.psiFull = 1.0 * std::atan2(spm.qyA + spm.qyC, spm.qxA + spm.qxC);
 
-    if (spm.centrality > cfg.cCentMax || spm.centrality < cfg.cCentMin)
+    if (spm.centrality > cfg.cCentMax || spm.centrality < cfg.cCentMin) {
       return;
+    }
 
     // Load correlations and SP resolution needed for Scalar Product and event plane methods.
     // Only load once!
     // If not loaded set to 1
 
-    if (cfg.cCCDBdir_QQ.value.empty() == false) {
+    if (!cfg.cCCDBdir_QQ.value.empty()) {
       if (!conf.clQQ) {
-        TList* hcorrList = ccdb->getForTimeStamp<TList>(cfg.cCCDBdir_QQ.value, bc.timestamp());
-        conf.hcorrQQ = reinterpret_cast<TProfile*>(hcorrList->FindObject("qAqCXY"));
-        conf.hcorrQQx = reinterpret_cast<TProfile*>(hcorrList->FindObject("qAqCX"));
-        conf.hcorrQQy = reinterpret_cast<TProfile*>(hcorrList->FindObject("qAqCY"));
+        auto* hcorrList = ccdb->getForTimeStamp<TList>(cfg.cCCDBdir_QQ.value, bc.timestamp());
+        conf.hcorrQQ = dynamic_cast<TProfile*>(hcorrList->FindObject("qAqCXY"));
+        conf.hcorrQQx = dynamic_cast<TProfile*>(hcorrList->FindObject("qAqCX"));
+        conf.hcorrQQy = dynamic_cast<TProfile*>(hcorrList->FindObject("qAqCY"));
         conf.clQQ = true;
       }
       spm.corrQQ = conf.hcorrQQ->GetBinContent(conf.hcorrQQ->FindBin(spm.centrality));
@@ -1766,19 +1885,20 @@ struct FlowSP {
     }
 
     double evPlaneRes = 1.;
-    if (cfg.cCCDBdir_SP.value.empty() == false) {
+    if (!cfg.cCCDBdir_SP.value.empty()) {
       if (!conf.clEvPlaneRes) {
         conf.hEvPlaneRes = ccdb->getForTimeStamp<TProfile>(cfg.cCCDBdir_SP.value, bc.timestamp());
         conf.clEvPlaneRes = true;
       }
       evPlaneRes = conf.hEvPlaneRes->GetBinContent(conf.hEvPlaneRes->FindBin(spm.centrality));
-      if (evPlaneRes < 0)
+      if (evPlaneRes < 0) {
         LOGF(fatal, "<Cos(PsiA-PsiC)> > 0 for centrality %.2f! Cannot determine resolution.. Change centrality ranges!!!", spm.centrality);
+      }
       evPlaneRes = std::sqrt(evPlaneRes);
     }
 
     spm.centWeight = 1.;
-    if (cfg.cCCDBdir_centrality.value.empty() == false) {
+    if (!cfg.cCCDBdir_centrality.value.empty()) {
       if (!conf.clCentrality) {
         conf.hCentrality = ccdb->getForTimeStamp<TH1D>(cfg.cCCDBdir_centrality.value, bc.timestamp());
         conf.clCentrality = true;
@@ -1799,8 +1919,9 @@ struct FlowSP {
 
       histos.fill(HIST("hPIDcounts"), trackPID, track.pt());
 
-      if (track.sign() == 0)
+      if (track.sign() == 0) {
         continue;
+      }
 
       histos.fill(HIST("hTrackCount"), trackSel_ZeroCharge);
 
@@ -1822,8 +1943,9 @@ struct FlowSP {
         }
       }
 
-      if (!trackSelected(track, field))
+      if (!trackSelected(track, field)) {
         continue;
+      }
 
       // constrain angle to 0 -> [0,0+2pi]
       auto phi = RecoDecay::constrainAngle(track.phi(), 0);
@@ -1913,7 +2035,7 @@ struct FlowSP {
 
   PROCESS_SWITCH(FlowSP, processDataPID, "Process analysis for non-derived data with PID", false);
 
-  void processMCReco(CC const& collision, aod::BCsWithTimestamps const&, TCs const& tracks, FilteredTCs const& filteredTracks, aod::McParticles const&)
+  void processMCReco(CC const& collision, aod::BCsWithTimestamps const&, TCs const& tracks, FilteredTCs const& filteredTracks, MCs const& McParts, aod::McCollisions const&)
   {
     auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
     int standardMagField = 99999;
@@ -1921,28 +2043,34 @@ struct FlowSP {
 
     spm.vz = collision.posZ();
     spm.centrality = collision.centFT0C();
-    if (cfg.cCentFT0Cvariant1)
+    if (cfg.cCentFT0Cvariant1) {
       spm.centrality = collision.centFT0CVariant1();
-    if (cfg.cCentFT0M)
+    }
+    if (cfg.cCentFT0M) {
       spm.centrality = collision.centFT0M();
-    if (cfg.cCentFV0A)
+    }
+    if (cfg.cCentFV0A) {
       spm.centrality = collision.centFV0A();
-    if (cfg.cCentNGlobal)
+    }
+    if (cfg.cCentNGlobal) {
       spm.centrality = collision.centNGlobal();
+    }
 
-    if (cfg.cFillQABefore)
+    if (cfg.cFillQABefore) {
       fillEventQA<kBefore>(collision, filteredTracks);
+    }
 
-    if (!eventSelected(collision, filteredTracks.size()))
+    if (!eventSelected(collision, filteredTracks.size())) {
       return;
+    }
 
-    if (spm.centrality > cfg.cCentMax || spm.centrality < cfg.cCentMin)
+    if (spm.centrality > cfg.cCentMax || spm.centrality < cfg.cCentMin) {
       return;
+    }
 
     histos.fill(HIST("hEventCount"), evSel_CentCuts);
 
     if (!collision.has_mcCollision()) {
-      LOGF(info, "No mccollision found for this collision");
       return;
     }
 
@@ -1951,15 +2079,17 @@ struct FlowSP {
     registry.fill(HIST("trackMCReco/hTrackSize_unFiltered"), tracks.size(), spm.centrality);
     registry.fill(HIST("trackMCReco/hTrackSize_Filtered"), filteredTracks.size(), spm.centrality);
 
-    for (const auto& track : filteredTracks) {
+    for (const auto& track : tracks) {
 
-      if (!track.has_mcParticle())
+      if (!track.has_mcParticle()) {
         continue;
+      }
 
-      auto mcParticle = track.mcParticle();
+      auto mcParticle = track.mcParticle_as<MCs>();
 
-      if (track.sign() == 0.0)
+      if (track.sign() == 0.0) {
         continue;
+      }
       histos.fill(HIST("hTrackCount"), trackSel_ZeroCharge);
 
       if (cfg.cFillWithMCParticle) {
@@ -1993,24 +2123,23 @@ struct FlowSP {
         }
       }
 
-      fillPrimaryHistos<kBefore, kInclusive>(mcParticle);
-      if (spm.charge == kPositive) {
-        fillPrimaryHistos<kBefore, kPositive>(mcParticle);
-      } else {
-        fillPrimaryHistos<kBefore, kNegative>(mcParticle);
+      if (!trackSelected(track, field)) {
+        continue;
       }
 
-      if (!trackSelected(track, field))
+      if (!mcParticle.isPhysicalPrimary()) {
         continue;
+      }
 
-      fillPrimaryHistos<kAfter, kInclusive>(mcParticle);
+      registry.fill(HIST("trackMCReco/incl/hPtMCPtTrack"), mcParticle.pt(), track.pt());
+      registry.fill(HIST("trackMCReco/incl/hEtaMCEtaTrack"), mcParticle.eta(), track.eta());
       if (spm.charge == kPositive) {
-        fillPrimaryHistos<kAfter, kPositive>(mcParticle);
-      } else {
-        fillPrimaryHistos<kAfter, kNegative>(mcParticle);
+        registry.fill(HIST("trackMCReco/pos/hPtMCPtTrack"), mcParticle.pt(), track.pt());
+        registry.fill(HIST("trackMCReco/pos/hEtaMCEtaTrack"), mcParticle.eta(), track.eta());
+      } else if (spm.charge == kNegative) {
+        registry.fill(HIST("trackMCReco/neg/hPtMCPtTrack"), mcParticle.pt(), track.pt());
+        registry.fill(HIST("trackMCReco/neg/hEtaMCEtaTrack"), mcParticle.eta(), track.eta());
       }
-      if (!mcParticle.isPhysicalPrimary())
-        continue;
 
       if (cfg.cFillWithMCParticle) {
         fillMCPtHistos<kAfter, kReco>(mcParticle, mcParticle.pdgCode());
@@ -2035,6 +2164,60 @@ struct FlowSP {
       }
 
     } // end of track loop
+
+    if (cfg.cDoGeneratedInReco) {
+
+      auto mcCollision = collision.mcCollision();
+      float vtxz = mcCollision.posZ();
+
+      // get McParticles which belong to mccollision
+      auto partSlice = McParts.sliceBy(partPerMcCollision, mcCollision.globalIndex());
+
+      for (const auto& particle : partSlice) {
+
+        if (!particle.isPhysicalPrimary()) {
+          continue;
+        }
+
+        auto pdgCode = particle.pdgCode();
+        auto pdgInfo = pdg->GetParticle(pdgCode);
+
+        if (std::abs(pdgInfo->Charge()) < 1) {
+          continue;
+        }
+
+        spm.charge = (pdgInfo->Charge() > 0) ? kPositive : kNegative;
+
+        int minVal = 100;
+        if (cfg.cFilterLeptons && std::abs(pdgCode) < minVal) {
+          continue;
+        }
+
+        fillMCPtHistos<kBefore, kGen>(particle, pdgCode);
+
+        registry.fill(HIST("trackMCGen/before/incl/phi_eta_vtxZ_gen"), particle.phi(), particle.eta(), vtxz);
+
+        if (spm.charge == kPositive) {
+          registry.fill(HIST("trackMCGen/before/pos/phi_eta_vtxZ_gen"), particle.phi(), particle.eta(), vtxz);
+        } else {
+          registry.fill(HIST("trackMCGen/before/neg/phi_eta_vtxZ_gen"), particle.phi(), particle.eta(), vtxz);
+        }
+
+        if (particle.eta() < -cfg.cTrackSelsEta || particle.eta() > cfg.cTrackSelsEta || particle.pt() < cfg.cTrackSelsPtmin || particle.pt() > cfg.cTrackSelsPtmax) {
+          continue;
+        }
+
+        fillMCPtHistos<kAfter, kGen>(particle, pdgCode);
+
+        registry.fill(HIST("trackMCGen/after/incl/phi_eta_vtxZ_gen"), particle.phi(), particle.eta(), vtxz);
+
+        if (spm.charge == kPositive) {
+          registry.fill(HIST("trackMCGen/after/pos/phi_eta_vtxZ_gen"), particle.phi(), particle.eta(), vtxz);
+        } else {
+          registry.fill(HIST("trackMCGen/after/neg/phi_eta_vtxZ_gen"), particle.phi(), particle.eta(), vtxz);
+        }
+      }
+    }
   }
   PROCESS_SWITCH(FlowSP, processMCReco, "Process analysis for MC reconstructed events", false);
 
@@ -2064,17 +2247,22 @@ struct FlowSP {
         auto filteredTrackSlice = filteredTracks.sliceBy(trackPerCollision, col.globalIndex());
 
         spm.centrality = col.centFT0C();
-        if (cfg.cCentFT0Cvariant1)
+        if (cfg.cCentFT0Cvariant1) {
           spm.centrality = col.centFT0CVariant1();
-        if (cfg.cCentFT0M)
+        }
+        if (cfg.cCentFT0M) {
           spm.centrality = col.centFT0M();
-        if (cfg.cCentFV0A)
+        }
+        if (cfg.cCentFV0A) {
           spm.centrality = col.centFV0A();
-        if (cfg.cCentNGlobal)
+        }
+        if (cfg.cCentNGlobal) {
           spm.centrality = col.centNGlobal();
+        }
 
-        if (cfg.cFillQABefore)
+        if (cfg.cFillQABefore) {
           fillEventQA<kBefore>(col, filteredTrackSlice);
+        }
 
         if (trackSlice.size() < 1) {
           colSelected = false;
@@ -2095,20 +2283,23 @@ struct FlowSP {
 
       } // leave reconstructed collision loop
 
-      if (!colSelected)
+      if (!colSelected) {
         continue;
+      }
 
       float vtxz = mcCollision.posZ();
 
       for (const auto& particle : partSlice) {
-        if (!particle.isPhysicalPrimary())
+        if (!particle.isPhysicalPrimary()) {
           continue;
+        }
 
         auto pdgCode = particle.pdgCode();
         auto pdgInfo = pdg->GetParticle(pdgCode);
 
-        if (std::abs(pdgInfo->Charge()) < 1)
+        if (std::abs(pdgInfo->Charge()) < 1) {
           continue;
+        }
 
         spm.charge = (pdgInfo->Charge() > 0) ? kPositive : kNegative;
 
@@ -2127,8 +2318,9 @@ struct FlowSP {
           registry.fill(HIST("trackMCGen/before/neg/phi_eta_vtxZ_gen"), particle.phi(), particle.eta(), vtxz);
         }
 
-        if (particle.eta() < -cfg.cTrackSelsEta || particle.eta() > cfg.cTrackSelsEta || particle.pt() < cfg.cTrackSelsPtmin || particle.pt() > cfg.cTrackSelsPtmax)
+        if (particle.eta() < -cfg.cTrackSelsEta || particle.eta() > cfg.cTrackSelsEta || particle.pt() < cfg.cTrackSelsPtmin || particle.pt() > cfg.cTrackSelsPtmax) {
           continue;
+        }
 
         fillMCPtHistos<kAfter, kGen>(particle, pdgCode);
 
@@ -2139,7 +2331,7 @@ struct FlowSP {
         } else {
           registry.fill(HIST("trackMCGen/after/neg/phi_eta_vtxZ_gen"), particle.phi(), particle.eta(), vtxz);
         }
-      }
+      } // end of parts
     }
   }
   PROCESS_SWITCH(FlowSP, processMCGen, "Process analysis for MC generated events", false);