Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mhemmer-cern merged 2 commits into from
Feb 21, 2026
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d921651
[PWGEM] PhotonMeson: Add second iteration of non lin for PCM
mhemmer-cern Feb 20, 2026
d3e553b
remove unsued variable
mhemmer-cern Feb 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 32 additions & 9 deletions PWGEM/PhotonMeson/Core/EMCPhotonCut.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
}
Expand All @@ -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]);
Expand All @@ -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;
}
Expand All @@ -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]);
Expand Down
71 changes: 40 additions & 31 deletions PWGEM/PhotonMeson/Core/EMNonLin.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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];
}
Loading
Loading