[PWGEM] PhotonMeson: Fix bug in EMNonLin

- Fixed a bug where for higher iterations of the NonLin correction the input for iteration i was not the input value * scale_(i-1) but the unscaled input value
- Changed the formluar of the NonLin from (1 + par0/x + par1/x^2) / (1 + par2/x) to (x + par0 +par1/x) / (x + par2) which is mathmatecially the same, however, it is a lot more stable since there is no 1/x^2 anymore.
- Optimized the `setParams` calling inside `nonLinProducer.cxx` where it now should only be called if the collision and hence the centrality changes. Before this it was called for each photon!

Author: mhemmer-cern

PWGEM/PhotonMeson/Core/EMNonLin.cxx

@@ -25,22 +25,27 @@ float EMNonLin::getCorrectionFactor(float x, const Context& ctx)
     return x;
   }
 
-  float val = x;
-  // safety measure
   int maxIter = std::min(ctx.nIter, MaxIter - 1);
 
+  float scale = 1.f; // cumulative scale
+  float refVal = x;  // reference value for computing next scale
+
   for (int i = 0; i <= maxIter; ++i) {
-    if (val == 0.f) {
+    if (refVal == 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);
-  }
+    const auto& p = ctx.params[i];
 
-  return val;
+    // scale function (x + a + b/x)/(x + c) which goes towards 1 for large x since x >> a,b,c -> x/x = 1
+    float iterScale =
+      (refVal + p.par0 + p.par1 / refVal) /
+      (refVal + p.par2);
+
+    scale *= iterScale; // total scale = product over itertaion scale
+    refVal = x * scale; // next iteration uses scaled original input
+  }
+  return scale;
 }
 
 const EMNonLin::NonLinParams* EMNonLin::resolveParams(PhotonType type, float cent)

PWGEM/PhotonMeson/TableProducer/nonLinProducer.cxx

@@ -109,32 +109,30 @@ struct NonLinProducer {
   template <o2::soa::is_table TClusters, o2::soa::is_iterator TCollisio>
   void runEMC(TClusters const& clusters, TCollisio& collision)
   {
-    float cent = getCentrality(collision);
 
     int32_t collIndex = collision.globalIndex();
+    float cent = getCentrality(collision);
+    emNonLinContextEMC.setParams(emNonLinEMC.resolveParams(o2::pwgem::nonlin::EMNonLin::PhotonType::kEMC, cent));
+
     for (const auto& cluster : clusters) {
-      float nonLinE = 0.f;
-      float nonLinPt = 0.f;
-      float nonLinFactor = 1.f;
 
       // check that we are at the correct collision
       if (cluster.emphotoneventId() != collIndex) {
         collIndex = cluster.emphotoneventId();
         collision.setCursor(collIndex);
         cent = getCentrality(collision);
+        emNonLinContextEMC.setParams(emNonLinEMC.resolveParams(o2::pwgem::nonlin::EMNonLin::PhotonType::kEMC, cent));
       }
 
-      emNonLinContextEMC.setParams(emNonLinEMC.resolveParams(o2::pwgem::nonlin::EMNonLin::PhotonType::kEMC, cent));
-
       // fill before non lin histograms
       historeg.fill(HIST("QA/EMC/EIn"), cluster.e());
       historeg.fill(HIST("QA/EMC/PtIn"), cluster.pt());
 
       // get NonLin factor from class dependent on the centrality
-      nonLinFactor = emNonLinEMC.getCorrectionFactor(cluster.e(), emNonLinContextEMC);
+      float nonLinFactor = emNonLinEMC.getCorrectionFactor(cluster.e(), emNonLinContextEMC);
 
-      nonLinE = nonLinFactor * cluster.e();
-      nonLinPt = nonLinFactor * cluster.pt();
+      float nonLinE = nonLinFactor * cluster.e();
+      float nonLinPt = nonLinFactor * cluster.pt();
 
       // fill after non lin histograms
       historeg.fill(HIST("QA/EMC/EOut"), nonLinE);
@@ -147,29 +145,27 @@ struct NonLinProducer {
   template <o2::soa::is_table TV0, o2::soa::is_iterator TCollisio>
   void runPCM(TV0 const& v0s, TCollisio& collision)
   {
-    float cent = getCentrality(collision);
 
     int32_t collIndex = collision.globalIndex();
+    float cent = getCentrality(collision);
+    emNonLinContextPCM.setParams(emNonLinPCM.resolveParams(o2::pwgem::nonlin::EMNonLin::PhotonType::kPCM, cent));
     for (const auto& v0 : v0s) {
-      float nonLinPt = 0.f;
-      float nonLinFactor = 1.f;
 
       // check that we are at the correct collision
       if (v0.emphotoneventId() != collIndex) {
         collIndex = v0.emphotoneventId();
         collision.setCursor(collIndex);
         cent = getCentrality(collision);
+        emNonLinContextPCM.setParams(emNonLinPCM.resolveParams(o2::pwgem::nonlin::EMNonLin::PhotonType::kPCM, cent));
       }
 
-      emNonLinContextPCM.setParams(emNonLinPCM.resolveParams(o2::pwgem::nonlin::EMNonLin::PhotonType::kPCM, cent));
-
       // fill before non lin histograms
       historeg.fill(HIST("QA/PCM/PtIn"), v0.pt());
 
       // get NonLin factor from class dependent on the centrality
-      nonLinFactor = emNonLinEMC.getCorrectionFactor(v0.pt(), emNonLinContextPCM);
+      float nonLinFactor = emNonLinEMC.getCorrectionFactor(v0.pt(), emNonLinContextPCM);
 
-      nonLinPt = nonLinFactor * v0.pt();
+      float nonLinPt = nonLinFactor * v0.pt();
 
       // fill after non lin histograms
       historeg.fill(HIST("QA/PCM/PtOut"), nonLinPt);