diff --git a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx
index 2e64bd6ebe2..3f019289988 100644
--- a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx
+++ b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx
@@ -14,69 +14,89 @@
 ///
 /// \author M. Hemmer, marvin.hemmer@cern.ch
 
-#include <numbers>
-#include <iterator>
-#include <array>
-#include <string>
-#include <map>
-#include <algorithm>
-#include <vector>
-#include <tuple>
-#include <utility>
-#include <cmath>
-
-#include "Math/Vector4D.h"
-#include "Math/Vector3D.h"
-#include "Math/LorentzRotation.h"
-#include "Math/Rotation3D.h"
-#include "Math/AxisAngle.h"
-
-#include "CCDB/BasicCCDBManager.h"
-#include "Framework/AnalysisTask.h"
-#include "Framework/ASoAHelpers.h"
-#include "Framework/HistogramRegistry.h"
-#include "Framework/runDataProcessing.h"
-
-#include "Common/Core/EventPlaneHelper.h"
-#include "Common/Core/RecoDecay.h"
-#include "Common/DataModel/Qvectors.h"
-#include "CommonConstants/MathConstants.h"
-
-#include "DetectorsBase/GeometryManager.h"
-#include "DataFormatsEMCAL/Constants.h"
-#include "EMCALBase/Geometry.h"
-#include "EMCALCalib/BadChannelMap.h"
-
-#include "PWGEM/Dilepton/Utils/EMTrackUtilities.h"
 #include "PWGEM/PhotonMeson/Core/EMCPhotonCut.h"
 #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h"
 #include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
-#include "PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h"
-#include "PWGEM/PhotonMeson/Utils/PairUtilities.h"
 #include "PWGEM/PhotonMeson/Utils/EventHistograms.h"
-#include "PWGEM/PhotonMeson/Utils/NMHistograms.h"
+#include "PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h"
+
+#include "Common/CCDB/TriggerAliases.h"
+#include "Common/Core/EventPlaneHelper.h"
+#include "Common/Core/RecoDecay.h"
+#include "Common/DataModel/Centrality.h"
+#include "Common/DataModel/EventSelection.h"
+
+#include <CCDB/BasicCCDBManager.h>
+#include <CommonConstants/MathConstants.h>
+#include <EMCALBase/Geometry.h>
+#include <EMCALBase/GeometryBase.h>
+#include <EMCALCalib/BadChannelMap.h>
+#include <Framework/ASoA.h>
+#include <Framework/ASoAHelpers.h>
+#include <Framework/AnalysisDataModel.h>
+#include <Framework/AnalysisHelpers.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/BinningPolicy.h>
+#include <Framework/Configurable.h>
+#include <Framework/Expressions.h>
+#include <Framework/GroupedCombinations.h>
+#include <Framework/HistogramRegistry.h>
+#include <Framework/HistogramSpec.h>
+#include <Framework/InitContext.h>
+#include <Framework/OutputObjHeader.h>
+#include <Framework/SliceCache.h>
+#include <Framework/runDataProcessing.h>
+
+#include <Math/GenVector/AxisAngle.h>
+#include <Math/GenVector/Rotation3D.h>
+#include <Math/Vector4D.h> // IWYU pragma: keep
+#include <TH1.h>
+#include <TString.h>
+
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <numbers>
+#include <string>
+#include <string_view>
+#include <tuple>
+#include <utility>
+#include <vector>
 
 using namespace o2;
 using namespace o2::aod;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
 using namespace o2::soa;
-using namespace o2::aod::pwgem::photonmeson::photonpair;
 using namespace o2::aod::pwgem::photon;
-using namespace o2::aod::pwgem::dilepton::utils;
-
-enum QvecEstimator { FT0M = 0,
-                     FT0A = 1,
-                     FT0C,
-                     TPCPos,
-                     TPCNeg,
-                     TPCTot };
-
-enum CentralityEstimator { None = 0,
-                           CFT0A = 1,
-                           CFT0C,
-                           CFT0M,
-                           NCentralityEstimators
+
+enum QvecEstimator {
+  FT0M = 0,
+  FT0A = 1,
+  FT0C,
+  TPCPos,
+  TPCNeg,
+  TPCTot
+};
+
+enum CentralityEstimator {
+  None = 0,
+  CFT0A = 1,
+  CFT0C,
+  CFT0M,
+  NCentralityEstimators
+};
+
+enum Harmonics {
+  kNone = 0,
+  kDirect = 1,
+  kElliptic = 2,
+  kTriangluar = 3,
+  kQuadrangular = 4,
+  kPentagonal = 5,
+  kHexagonal = 6,
+  kHeptagonal = 7,
+  kOctagonal = 8
 };
 
 struct TaskPi0FlowEMC {
@@ -98,6 +118,8 @@ struct TaskPi0FlowEMC {
   Configurable<bool> cfgDoM02{"cfgDoM02", false, "Flag to enable flow vs M02 for single photons"};
   Configurable<bool> cfgDoReverseScaling{"cfgDoReverseScaling", false, "Flag to reverse the scaling that is possibly applied during NonLin"};
   Configurable<bool> cfgDoPlaneQA{"cfgDoPlaneQA", false, "Flag to enable QA plots comparing in and out of plane"};
+  Configurable<float> cfgMaxQVector{"cfgMaxQVector", 20.f, "Maximum allowed absolute QVector value."};
+  Configurable<float> cfgMaxAsymmetry{"cfgMaxAsymmetry", 0.1f, "Maximum allowed asymmetry for photon pairs used in calibration."};
 
   // configurable axis
   ConfigurableAxis thnConfigAxisInvMass{"thnConfigAxisInvMass", {400, 0.0, 0.8}, ""};
@@ -205,6 +227,9 @@ struct TaskPi0FlowEMC {
   std::vector<int8_t> lookupTable1D;
   float epsilon = 1.e-8;
 
+  // static constexpr
+  static constexpr int64_t NMinPhotonRotBkg = 3;
+
   // To access the 1D array
   inline int getIndex(int iEta, int iPhi)
   {
@@ -265,7 +290,7 @@ struct TaskPi0FlowEMC {
 
   void init(InitContext&)
   {
-    if (harmonic != 2 && harmonic != 3) {
+    if (harmonic != kElliptic && harmonic != kTriangluar) {
       LOG(info) << "Harmonic was set to " << harmonic << " but can only be 2 or 3!";
     }
 
@@ -492,65 +517,65 @@ struct TaskPi0FlowEMC {
 
     switch (detector) {
       case QvecEstimator::FT0M:
-        if (harmonic == 2) {
+        if (harmonic == kElliptic) {
           xQVec = collision.q2xft0m();
           yQVec = collision.q2yft0m();
-        } else if (harmonic == 3) {
+        } else if (harmonic == kTriangluar) {
           xQVec = collision.q3xft0m();
           yQVec = collision.q3yft0m();
         }
         break;
       case QvecEstimator::FT0A:
-        if (harmonic == 2) {
+        if (harmonic == kElliptic) {
           xQVec = collision.q2xft0a();
           yQVec = collision.q2yft0a();
-        } else if (harmonic == 3) {
+        } else if (harmonic == kTriangluar) {
           xQVec = collision.q3xft0a();
           yQVec = collision.q3yft0a();
         }
         break;
       case QvecEstimator::FT0C:
-        if (harmonic == 2) {
+        if (harmonic == kElliptic) {
           xQVec = collision.q2xft0c();
           yQVec = collision.q2yft0c();
-        } else if (harmonic == 3) {
+        } else if (harmonic == kTriangluar) {
           xQVec = collision.q3xft0c();
           yQVec = collision.q3yft0c();
         }
         break;
       case QvecEstimator::TPCPos:
-        if (harmonic == 2) {
+        if (harmonic == kElliptic) {
           xQVec = collision.q2xbpos();
           yQVec = collision.q2ybpos();
-        } else if (harmonic == 3) {
+        } else if (harmonic == kTriangluar) {
           xQVec = collision.q3xbpos();
           yQVec = collision.q3ybpos();
         }
         break;
       case QvecEstimator::TPCNeg:
-        if (harmonic == 2) {
+        if (harmonic == kElliptic) {
           xQVec = collision.q2xbneg();
           yQVec = collision.q2ybneg();
-        } else if (harmonic == 3) {
+        } else if (harmonic == kTriangluar) {
           xQVec = collision.q3xbneg();
           yQVec = collision.q3ybneg();
         }
         break;
       case QvecEstimator::TPCTot:
-        if (harmonic == 2) {
+        if (harmonic == kElliptic) {
           xQVec = collision.q2xbtot();
           yQVec = collision.q2ybtot();
-        } else if (harmonic == 3) {
+        } else if (harmonic == kTriangluar) {
           xQVec = collision.q3xbtot();
           yQVec = collision.q3ybtot();
         }
         break;
       default:
         LOG(warning) << "Q vector estimator not valid. Falling back to FT0M";
-        if (harmonic == 2) {
+        if (harmonic == kElliptic) {
           xQVec = collision.q2xft0m();
           yQVec = collision.q2yft0m();
-        } else if (harmonic == 3) {
+        } else if (harmonic == kTriangluar) {
           xQVec = collision.q3xft0m();
           yQVec = collision.q3yft0m();
         }
@@ -565,7 +590,7 @@ struct TaskPi0FlowEMC {
   {
     bool isgood = true;
     for (const auto& QVec : QVecs) {
-      if (std::fabs(QVec) > 20.f) {
+      if (std::fabs(QVec) > cfgMaxQVector) {
         isgood = false;
         break;
       }
@@ -660,7 +685,7 @@ struct TaskPi0FlowEMC {
   void rotationBackground(const ROOT::Math::PtEtaPhiMVector& meson, ROOT::Math::PtEtaPhiMVector photon1, ROOT::Math::PtEtaPhiMVector photon2, TPhotons const& photons_coll, unsigned int ig1, unsigned int ig2, CollsWithQvecs::iterator const& collision)
   {
     // if less than 3 clusters are present skip event since we need at least 3 clusters
-    if (photons_coll.size() < 3) {
+    if (photons_coll.size() < NMinPhotonRotBkg) {
       return;
     }
 
@@ -759,7 +784,7 @@ struct TaskPi0FlowEMC {
   void rotationBackgroundCalib(const ROOT::Math::PtEtaPhiMVector& meson, ROOT::Math::PtEtaPhiMVector photon1, ROOT::Math::PtEtaPhiMVector photon2, TPhotons const& photons_coll, unsigned int ig1, unsigned int ig2, CollsWithQvecs::iterator const& collision)
   {
     // if less than 3 clusters are present skip event since we need at least 3 clusters
-    if (photons_coll.size() < 3) {
+    if (photons_coll.size() < NMinPhotonRotBkg) {
       return;
     }
     float cent = getCentrality(collision);
@@ -794,7 +819,7 @@ struct TaskPi0FlowEMC {
       }
       ROOT::Math::PtEtaPhiMVector photon3(photon.pt(), photon.eta(), photon.phi(), 0.);
       if (iCellIDPhoton1 >= 0) {
-        if (std::fabs((photon1.E() - photon3.E()) / (photon1.E() + photon3.E()) < 0.1)) { // only use symmetric decays
+        if (std::fabs((photon1.E() - photon3.E()) / (photon1.E() + photon3.E()) < cfgMaxAsymmetry)) { // only use symmetric decays
           ROOT::Math::PtEtaPhiMVector mother1 = photon1 + photon3;
           float openingAngle1 = std::acos(photon1.Vect().Dot(photon3.Vect()) / (photon1.P() * photon3.P()));
 
@@ -812,7 +837,7 @@ struct TaskPi0FlowEMC {
         }
       }
       if (iCellIDPhoton2 >= 0) {
-        if (std::fabs((photon2.E() - photon3.E()) / (photon2.E() + photon3.E()) < 0.1)) { // only use symmetric decays
+        if (std::fabs((photon2.E() - photon3.E()) / (photon2.E() + photon3.E()) < cfgMaxAsymmetry)) { // only use symmetric decays
           ROOT::Math::PtEtaPhiMVector mother2 = photon2 + photon3;
           float openingAngle2 = std::acos(photon2.Vect().Dot(photon3.Vect()) / (photon2.P() * photon3.P()));
 
@@ -1138,7 +1163,7 @@ struct TaskPi0FlowEMC {
     float yQVecBNeg = -999.f;
     float xQVecBTot = -999.f;
     float yQVecBTot = -999.f;
-    if (harmonic == 2) {
+    if (harmonic == kElliptic) {
       xQVecFT0a = collision.q2xft0a();
       yQVecFT0a = collision.q2yft0a();
       xQVecFT0c = collision.q2xft0c();
@@ -1151,7 +1176,7 @@ struct TaskPi0FlowEMC {
       yQVecBNeg = collision.q2ybneg();
       xQVecBTot = collision.q2xbtot();
       yQVecBTot = collision.q2ybtot();
-    } else if (harmonic == 3) {
+    } else if (harmonic == kTriangluar) {
       xQVecFT0a = collision.q3xft0a();
       yQVecFT0a = collision.q3yft0a();
       xQVecFT0c = collision.q3xft0c();
@@ -1203,7 +1228,7 @@ struct TaskPi0FlowEMC {
       registry.fill(HIST("epReso/hEpResoFT0mTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0m, epBNegs)));
       registry.fill(HIST("epReso/hEpResoFT0mTPCtot"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0m, epBTots)));
       registry.fill(HIST("epReso/hEpResoTPCposTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epBPoss, epBNegs)));
-      for (int n = 1; n <= 8; n++) {
+      for (int n = 1; n <= kOctagonal; n++) {
         registry.fill(HIST("epReso/hEpCosCoefficientsFT0c"), centrality, n, std::cos(n * epFT0c));
         registry.fill(HIST("epReso/hEpSinCoefficientsFT0c"), centrality, n, std::sin(n * epFT0c));
         registry.fill(HIST("epReso/hEpCosCoefficientsFT0a"), centrality, n, std::cos(n * epFT0a));
@@ -1323,7 +1348,7 @@ struct TaskPi0FlowEMC {
           registry.fill(HIST("hClusterCuts"), 5);
           continue;
         }
-        if (std::fabs((v1.E() - v2.E()) / (v1.E() + v2.E()) < 0.1)) { // only use symmetric decays
+        if (std::fabs((v1.E() - v2.E()) / (v1.E() + v2.E()) < cfgMaxAsymmetry)) { // only use symmetric decays
           registry.fill(HIST("hClusterCuts"), 6);
           registry.fill(HIST("hSparseCalibSE"), vMeson.M(), vMeson.E() / 2., getCentrality(collision));
         }