[PWGEM] PhotonMeson: Add second iteration of NonLin for PCM (#15073)

Author: mhemmer-cern

PWGEM/PhotonMeson/Core/EMCPhotonCut.h

@@ -41,6 +41,20 @@ static constexpr bool HasPrimaries = !std::is_same_v<TPrimaries, std::nullptr_t>
 template <typename TSecondaries>
 static constexpr bool HasSecondaries = !std::is_same_v<TSecondaries, std::nullptr_t>;
 
+template <typename T>
+concept IsNonLinIterator = o2::soa::is_iterator<T> && requires(T t) {
+  // Check that the *elements* of the container have the required methods:
+  { t.corrE() } -> std::same_as<float>;
+  { t.corrPt() } -> std::same_as<float>;
+};
+
+template <typename T>
+concept IsNonLinContainer = o2::soa::is_table<T> && requires(T t) {
+  // Check that the *elements* of the container have the required methods:
+  { t.begin().corrE() } -> std::same_as<float>;
+  { t.begin().corrPt() } -> std::same_as<float>;
+};
+
 template <typename T>
 concept IsTrackIterator = o2::soa::is_iterator<T> && requires(T t) {
   // Check that the *elements* of the container have the required methods:
@@ -114,7 +128,7 @@ class EMCPhotonCut : public TNamed
 
   static const char* mCutNames[static_cast<int>(EMCPhotonCuts::kNCuts)];
 
-  constexpr auto getClusterId(o2::soa::is_iterator auto const& t) const
+  static constexpr auto getClusterId(o2::soa::is_iterator auto const& t)
   {
     if constexpr (requires { t.emEmcClusterId(); }) {
       return t.emEmcClusterId();
@@ -384,7 +398,11 @@ class EMCPhotonCut : public TNamed
         return cluster.definition() == mDefinition;
 
       case EMCPhotonCuts::kEnergy:
-        return cluster.e() > mMinE;
+        if constexpr (IsNonLinIterator<std::decay_t<decltype(cluster)>>) {
+          return cluster.corrE() > mMinE;
+        } else {
+          return cluster.e() > mMinE;
+        }
 
       case EMCPhotonCuts::kNCell:
         return cluster.nCells() >= mMinNCell;
@@ -478,7 +496,11 @@ class EMCPhotonCut : public TNamed
         return cluster.definition() == mDefinition;
 
       case EMCPhotonCuts::kEnergy:
-        return cluster.e() > mMinE;
+        if constexpr (IsNonLinIterator<Cluster>) {
+          return cluster.corrE() > mMinE;
+        } else {
+          return cluster.e() > mMinE;
+        }
 
       case EMCPhotonCuts::kNCell:
         return cluster.nCells() >= mMinNCell;
@@ -496,7 +518,9 @@ class EMCPhotonCut : public TNamed
             auto dPhi = std::fabs(emcmatchedtrack.deltaPhi());
             auto trackpt = emcmatchedtrack.trackPt();
             auto trackp = emcmatchedtrack.trackP();
-            bool result = (dEta > GetTrackMatchingEta(trackpt)) || (dPhi > GetTrackMatchingPhi(trackpt)) || (cluster.e() / trackp >= mMinEoverP);
+            bool result = (dEta > GetTrackMatchingEta(trackpt)) ||
+                          (dPhi > GetTrackMatchingPhi(trackpt)) ||
+                          (cluster.e() / trackp >= mMinEoverP);
             if (!result) {
               return false;
             }
@@ -506,7 +530,7 @@ class EMCPhotonCut : public TNamed
           auto dPhis = cluster.deltaPhi();   // std:vector<float>
           auto trackspt = cluster.trackpt(); // std:vector<float>
           auto tracksp = cluster.trackp();   // std:vector<float>
-          int ntrack = tracksp.size();
+          int ntrack = trackspt.size();
           for (int itr = 0; itr < ntrack; itr++) {
             float dEta = std::fabs(dEtas[itr]);
             float dPhi = std::fabs(dPhis[itr]);
@@ -526,8 +550,8 @@ class EMCPhotonCut : public TNamed
             auto dEta = std::fabs(emcmatchedtrack.deltaEta());
             auto dPhi = std::fabs(emcmatchedtrack.deltaPhi());
             auto trackpt = emcmatchedtrack.trackPt();
-            auto trackp = emcmatchedtrack.trackP();
-            bool result = (dEta > GetSecTrackMatchingEta(trackpt)) || (dPhi > GetSecTrackMatchingPhi(trackpt)) || (cluster.e() / trackp >= mMinEoverP);
+            bool result = (dEta > GetSecTrackMatchingEta(trackpt)) ||
+                          (dPhi > GetSecTrackMatchingPhi(trackpt));
             if (!result) {
               return false;
             }
@@ -536,8 +560,7 @@ class EMCPhotonCut : public TNamed
           auto dEtas = cluster.deltaEtaSec();   // std:vector<float>
           auto dPhis = cluster.deltaPhiSec();   // std:vector<float>
           auto trackspt = cluster.trackptSec(); // std:vector<float>
-          auto tracksp = cluster.trackpSec();   // std:vector<float>
-          int ntrack = tracksp.size();
+          int ntrack = trackspt.size();
           for (int itr = 0; itr < ntrack; itr++) {
             float dEta = std::fabs(dEtas[itr]);
             float dPhi = std::fabs(dPhis[itr]);

PWGEM/PhotonMeson/Core/EMNonLin.cxx

@@ -15,43 +15,52 @@
 
 #include "EMNonLin.h"
 
-#include <cmath>
+#include <algorithm>
 
 using namespace o2::pwgem::nonlin;
 
-float EMNonLin::getCorrectionFactor(float inputCalibValue, PhotonType photonType, float cent)
+float EMNonLin::getCorrectionFactor(float x, const Context& ctx)
 {
-  float param0 = 0, param1 = 0, param2 = 0, val = 1.f;
-  switch (photonType) {
-    case PhotonType::kEMC:
-      if (cent >= 30 && cent <= 40) {
-        param0 = -5.33426e-01f;
-        param1 = 1.40144e-02f;
-        param2 = -5.24434e-01f;
-      } else {
-        param0 = 0.f;
-        param1 = 0.f;
-        param2 = 0.f;
-      }
-      break;
-    case PhotonType::kPCM:
-      if (cent >= 0 && cent <= 100) {
-        param0 = 10.7203f;
-        param1 = 0.0383968f;
-        param2 = 10.6025f;
-      } else {
-        param0 = 0.f;
-        param1 = 0.f;
-        param2 = 0.f;
-      }
-      break;
-    case PhotonType::kPHOS:
-      param0 = 0.f;
-      param1 = 0.f;
-      param2 = 0.f;
+  if (!ctx.params || x == 0.f) [[unlikely]] {
+    return x;
+  }
+
+  float val = x;
+  // safety measure
+  int maxIter = std::min(ctx.nIter, MaxIter - 1);
+
+  for (int i = 0; i <= maxIter; ++i) {
+    if (val == 0.f) {
       break;
+    }
+
+    float inv = 1.f / val;
+    val *= (1.f + ctx.params[i].par0 * inv +
+            ctx.params[i].par1 * inv * inv) /
+           (1.f + ctx.params[i].par2 * inv);
   }
 
-  val = (1.f + param0 / inputCalibValue + param1 / std::pow(inputCalibValue, 2.f)) / (1.f + param2 / inputCalibValue);
   return val;
 }
+
+const EMNonLin::NonLinParams* EMNonLin::resolveParams(PhotonType type, float cent)
+{
+  int centBin = static_cast<int>(cent / 10.f);
+  if (centBin < 0)
+    centBin = 0;
+  if (centBin >= CentBins)
+    centBin = CentBins - 1;
+
+  return &kNonLinTable[static_cast<int>(type)][centBin][0];
+}
+
+const EMNonLin::NonLinParams* EMNonLin::resolveParamsMC(PhotonType type, float cent)
+{
+  int centBin = static_cast<int>(cent / 10.f);
+  if (centBin < 0)
+    centBin = 0;
+  if (centBin >= CentBins)
+    centBin = CentBins - 1;
+
+  return &kNonLinTableMC[static_cast<int>(type)][centBin][0];
+}