diff --git a/PWGEM/PhotonMeson/Core/EMNonLin.cxx b/PWGEM/PhotonMeson/Core/EMNonLin.cxx index bb5102dfc2a..571d7efbae8 100644 --- a/PWGEM/PhotonMeson/Core/EMNonLin.cxx +++ b/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) diff --git a/PWGEM/PhotonMeson/TableProducer/nonLinProducer.cxx b/PWGEM/PhotonMeson/TableProducer/nonLinProducer.cxx index 9d3847ffb8a..74ac6e4ed5e 100644 --- a/PWGEM/PhotonMeson/TableProducer/nonLinProducer.cxx +++ b/PWGEM/PhotonMeson/TableProducer/nonLinProducer.cxx @@ -109,32 +109,30 @@ struct NonLinProducer { template 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 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);