Skip to content

[PWGLF] Add new process function for MC spin correlation, fix systematic selection for f1-p correlation #15230

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
skundu692 merged 9 commits into from
Mar 2, 2026
+726 −89
Merged

[PWGLF] Add new process function for MC spin correlation, fix systematic selection for f1-p correlation #15230

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
5077cd1
Add process function for systematic
skundu692 Jan 30, 2026
95bdcd0
Fix clang issue
skundu692 Feb 1, 2026
cade617
Systematic variation f1p
skundu692 Feb 24, 2026
a9fcaa7
Old process function
skundu692 Feb 25, 2026
edfdf66
Fill QA histogram for trigger events
skundu692 Feb 28, 2026
dcbdb80
Add Event loss in phi meson task and add MC process function for spin…
skundu692 Mar 2, 2026
07f1c68
Add new MC process function for mixing
skundu692 Mar 2, 2026
9cedb5a
MC process function
skundu692 Mar 2, 2026
b5d2e4f
Fix clang error
skundu692 Mar 2, 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
62 changes: 52 additions & 10 deletions PWGLF/TableProducer/Resonances/f1protonreducedtable.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,7 @@ struct f1protonreducedtable {
Configurable<int> trackSphMin{"trackSphMin", 10, "Number of tracks for Spherocity Calculation"};

// Configs for track PID
Configurable<bool> Confglobaltrackcheck{"Confglobaltrackcheck", true, "Global track check"};
Configurable<bool> cfgSkimmedProcessing{"cfgSkimmedProcessing", true, "Analysed skimmed events"};
Configurable<bool> ConfUseManualPIDproton{"ConfUseManualPIDproton", true, "True: use home-made PID solution for proton "};
Configurable<bool> ConfUseManualPIDkaon{"ConfUseManualPIDkaon", true, "True: use home-made PID solution for kaon "};
Expand Down Expand Up @@ -161,16 +162,16 @@ struct f1protonreducedtable {
{"hInvMassf1Like", "hInvMassf1Like", {HistType::kTH2F, {{400, 1.1f, 1.9f}, {100, 0.0f, 10.0f}}}},
{"hInvMassf1kstar", "hInvMassf1kstar", {HistType::kTH3F, {{400, 1.1f, 1.9f}, {100, 0.0f, 10.0f}, {8, 0.0f, 0.8f}}}},
{"hkstarDist", "hkstarDist", {HistType::kTH1F, {{300, 0.0f, 3.0f}}}},
{"hDCAxy", "hDCAxy", {HistType::kTH1F, {{100, -5.0f, 5.0f}}}},
{"hDCAz", "hDCAz", {HistType::kTH1F, {{100, -5.0f, 5.0f}}}},
{"hDCAxy", "hDCAxy", {HistType::kTH3F, {{100, -0.05f, 0.05f}, {5, -2.5, 2.5}, {40, 0.0, 4.0}}}},
{"hDCAz", "hDCAz", {HistType::kTH3F, {{100, -0.05f, 0.05f}, {2, 0, 2}, {40, 0.0, 4.0}}}},
{"hPhi", "hPhi", {HistType::kTH1F, {{1400, -7.0f, 7.0f}}}},
{"hPhiSphero", "hPhiSphero", {HistType::kTH1F, {{1400, -7.0f, 7.0f}}}},
{"hEta", "hEta", {HistType::kTH1F, {{20, -1.0f, 1.0f}}}},
{"hNsigmaPtpionTPC", "hNsigmaPtpionTPC", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
{"hNsigmaPtpionTOF", "hNsigmaPtpionTOF", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
{"hNsigmaPtkaonTPC", "hNsigmaPtkaonTPC", {HistType::kTH3F, {{200, -10.0f, 10.0f}, {200, -20.0f, 20.0f}, {100, 0.0f, 10.0f}}}},
{"hNsigmaPtkaonTOF", "hNsigmaPtkaonTOF", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
{"hNsigmaPtprotonTPC", "hNsigmaPtprotonTPC", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
{"hNsigmaPtprotonTPC", "hNsigmaPtprotonTPC", {HistType::kTH3F, {{100, -5.0f, 5.0f}, {100, -5.0f, 5.0f}, {100, 0.0f, 10.0f}}}},
{"hNsigmaPtprotonTOF", "hNsigmaPtprotonTOF", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
{"hInvMassk0", "hInvMassk0", {HistType::kTH2F, {{200, 0.4f, 0.6f}, {100, 0.0f, 10.0f}}}},
},
Expand Down Expand Up @@ -265,7 +266,7 @@ struct f1protonreducedtable {
template <typename T>
bool selectionGlobalTrack(const T& candidate)
{
if (!(candidate.isGlobalTrack() && candidate.isPVContributor())) {
if (Confglobaltrackcheck && !(candidate.isGlobalTrack() && candidate.isPVContributor())) {
return false;
}
return true;
Expand Down Expand Up @@ -348,6 +349,32 @@ struct f1protonreducedtable {
return false;
}

inline bool passProtonPID(float nsigmaTPC, float nsigmaTOF, float TOFHit, ROOT::Math::PtEtaPhiMVector proton)
{
// pidMode:
// 0 = default: p < thr -> |TPC| < 2.5 ; p >= thr -> TOF mandatory AND circular(TPC,TOF) < 2.0
// 1 = syst-1: p < thr -> |TPC| < 2.0 ; p >= thr -> TOF mandatory AND circular(TPC,TOF) < 2.0
// 2 = syst-2: p < thr -> |TPC| < 2.5 ; p >= thr -> TOF mandatory AND circular(TPC,TOF) < 2.5

if (proton.Pt() > 4.0 || proton.Pt() < 0.15) {
return false;
}

if (proton.P() < 0.7) {
return (std::abs(nsigmaTPC) < 2.5);
}

// above threshold: TOF must exist
if (TOFHit != 1) {
return false;
}

const float nsTPC = nsigmaTPC;
const float nsTOF = nsigmaTOF;
const float comb = std::sqrt(nsTPC * nsTPC + nsTOF * nsTOF);
return (comb < 2.5);
}

template <typename Collision, typename V0>
bool SelectionV0(Collision const& collision, V0 const& candidate)
{
Expand Down Expand Up @@ -682,8 +709,6 @@ struct f1protonreducedtable {
if (!selectionGlobalTrack(track)) {
continue;
}
qaRegistry.fill(HIST("hDCAxy"), track.dcaXY());
qaRegistry.fill(HIST("hDCAz"), track.dcaZ());
qaRegistry.fill(HIST("hEta"), track.eta());
qaRegistry.fill(HIST("hPhi"), track.phi());
double nTPCSigmaP[3]{track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()};
Expand Down Expand Up @@ -778,17 +803,15 @@ struct f1protonreducedtable {
ProtonIndex.push_back(track.globalIndex());
ProtonCharge.push_back(track.sign());

ProtonDcaxy.push_back(std::abs(track.dcaXY()));
ProtonDcaz.push_back(std::abs(track.dcaZ()));
ProtonDcaxy.push_back(track.dcaXY());
ProtonDcaz.push_back(track.dcaZ());
ProtonTPCNcls.push_back(std::abs(track.tpcNClsFound()));
ProtonTPCNcrs.push_back(std::abs(track.tpcNClsCrossedRows()));

if (track.sign() > 0) {
qaRegistry.fill(HIST("hNsigmaPtprotonTPC"), nTPCSigmaP[2], track.pt());
ProtonTPCNsigma.push_back(nTPCSigmaP[2]);
}
if (track.sign() < 0) {
qaRegistry.fill(HIST("hNsigmaPtprotonTPC"), nTPCSigmaN[2], track.pt());
ProtonTPCNsigma.push_back(nTPCSigmaN[2]);
}
if (track.hasTOF()) {
Expand Down Expand Up @@ -961,6 +984,25 @@ struct f1protonreducedtable {
}
}
qaRegistry.fill(HIST("hEventstat"), 0.5);
if (keepEventF1Proton) {
for (auto iproton = protons.begin(); iproton != protons.end(); ++iproton) {
auto i6 = std::distance(protons.begin(), iproton);
ProtonVectorDummy2 = protons.at(i6);
if (std::abs(ProtonDcaxy.at(i6)) < 0.05 && std::abs(ProtonDcaz.at(i6)) < 0.05) {
if (ProtonTOFHit.at(i6) && ProtonVectorDummy2.P() > 0.7) {
qaRegistry.fill(HIST("hNsigmaPtprotonTPC"), ProtonTPCNsigma.at(i6), ProtonTOFNsigma.at(i6), ProtonVectorDummy2.Pt());
}
if (ProtonVectorDummy2.P() < 0.7) {
qaRegistry.fill(HIST("hNsigmaPtprotonTPC"), ProtonTPCNsigma.at(i6), 4.999, ProtonVectorDummy2.Pt());
}
}
if (passProtonPID(ProtonTPCNsigma.at(i6), ProtonTOFNsigma.at(i6), ProtonTOFHit.at(i6), ProtonVectorDummy2)) {
qaRegistry.fill(HIST("hDCAxy"), ProtonDcaxy.at(i6), ProtonCharge.at(i6), ProtonVectorDummy2.Pt());
qaRegistry.fill(HIST("hDCAz"), ProtonDcaz.at(i6), ProtonCharge.at(i6), ProtonVectorDummy2.Pt());
}
}
}

if (numberF1 > 0 && (f1resonance.size() == f1signal.size()) && (f1resonance.size() == f1kaonkshortmass.size()) && (f1resonance.size() == f1resonanced1.size()) && (f1resonance.size() == f1resonanced2.size()) && (f1resonance.size() == f1resonanced3.size())) {
qaRegistry.fill(HIST("hEventstat"), 1.5);
if (keepEventF1Proton) {
Expand Down
156 changes: 84 additions & 72 deletions PWGLF/Tasks/Resonances/f1protoncorrelation.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
#include <random>
#include <sstream>
#include <string>
#include <unordered_set>
#include <vector>

using namespace o2;
Expand Down Expand Up @@ -140,86 +141,94 @@ struct f1protoncorrelation {
// -------------------------
void buildSystematicCuts()
{
// choices from your table/picture
const std::array<float, 2> optDcaxy{0.01f, 0.03f};
const std::array<float, 2> optDcaz{0.01f, 0.03f};
const std::array<int, 2> optNcrs{90, 100};
const std::array<int, 2> optNcls{90, 100};

const std::array<float, 2> optCPA{0.99f, 0.995f};
const std::array<float, 2> optRad{0.8f, 1.0f};
const std::array<float, 2> optDcaDD{0.9f, 0.8f};
const std::array<float, 2> optDcaV0{0.2f, 0.15f};
const std::array<float, 2> optLife{16.f, 18.f};
const std::array<float, 2> optDcaD1{0.06f, 0.08f};
const std::array<float, 2> optDcaD2{0.06f, 0.08f};

// bit layout:
// primary: 4 bits (0..3)
// v0: 7 bits (4..10)
// pid: 3 bits (11..13)
// total: 14 bits -> 16384 combos
auto buildFromMask = [&](uint32_t m) -> SysCuts {
// 3 options per cut: index 0 = DEFAULT, index 1/2 = variations
// Fill these with the exact values you want (I used your def as index 0 + your old options as 1/2)
const std::array<float, 3> optDcaxy{0.05f, 0.01f, 0.03f};
const std::array<float, 3> optDcaz{0.05f, 0.01f, 0.03f};
const std::array<int, 3> optNcrs{80, 90, 100};
const std::array<int, 3> optNcls{80, 90, 100};

const std::array<float, 3> optCPA{0.985f, 0.99f, 0.995f};
const std::array<float, 3> optRad{0.50f, 0.80f, 1.00f};
const std::array<float, 3> optDcaDD{1.00f, 0.90f, 0.80f};
const std::array<float, 3> optDcaV0{0.30f, 0.20f, 0.15f};
const std::array<float, 3> optLife{20.f, 16.f, 18.f};
const std::array<float, 3> optDcaD1{0.05f, 0.06f, 0.08f};
const std::array<float, 3> optDcaD2{0.05f, 0.06f, 0.08f};

// PID modes: also allow default (0) to appear in random variations
const std::array<int, 3> optPidPi{0, 1, 2};
const std::array<int, 3> optPidK{0, 1, 2};
const std::array<int, 3> optPidP{0, 1, 2};

// Helper: build SysCuts from chosen indices (0..2)
auto buildFromIdx = [&](int i0, int i1, int i2, int i3,
int i4, int i5, int i6, int i7, int i8, int i9, int i10,
int i11, int i12, int i13) -> SysCuts {
SysCuts c{};
c.maxDcaxy = optDcaxy[(m >> 0) & 1u];
c.maxDcaz = optDcaz[(m >> 1) & 1u];
c.minTPCCrossedRows = optNcrs[(m >> 2) & 1u];
c.minTPCClusters = optNcls[(m >> 3) & 1u];

c.minCPA = optCPA[(m >> 4) & 1u];
c.minRadius = optRad[(m >> 5) & 1u];
c.maxDcaDaughters = optDcaDD[(m >> 6) & 1u];
c.maxDcaV0 = optDcaV0[(m >> 7) & 1u];
c.maxLifetime = optLife[(m >> 8) & 1u];
c.minDcaD1 = optDcaD1[(m >> 9) & 1u];
c.minDcaD2 = optDcaD2[(m >> 10) & 1u];

c.pidPi = (m >> 11) & 1u;
c.pidK = (m >> 12) & 1u;
c.pidP = (m >> 13) & 1u;
c.maxDcaxy = optDcaxy[i0];
c.maxDcaz = optDcaz[i1];
c.minTPCCrossedRows = optNcrs[i2];
c.minTPCClusters = optNcls[i3];

c.minCPA = optCPA[i4];
c.minRadius = optRad[i5];
c.maxDcaDaughters = optDcaDD[i6];
c.maxDcaV0 = optDcaV0[i7];
c.maxLifetime = optLife[i8];
c.minDcaD1 = optDcaD1[i9];
c.minDcaD2 = optDcaD2[i10];

c.pidPi = optPidPi[i11];
c.pidK = optPidK[i12];
c.pidP = optPidP[i13];
return c;
};

// DEFAULT (sysId=0): set to your baseline
SysCuts def{};
def.maxDcaxy = 0.05f;
def.maxDcaz = 0.05f;
def.minTPCCrossedRows = 80;
def.minTPCClusters = 80;

def.minCPA = 0.985f;
def.minRadius = 0.5f;
def.maxDcaDaughters = 1.0f;
def.maxDcaV0 = 0.3f;
def.maxLifetime = 20.f;
def.minDcaD1 = 0.05f;
def.minDcaD2 = 0.05f;

def.pidPi = 0;
def.pidK = 0;
def.pidP = 0;

std::vector<uint32_t> masks;
masks.reserve(16384);
for (uint32_t m = 0; m < 16384; ++m) {
masks.push_back(m);
}
// sysId=0 must be strict default (all indices = 0)
SysCuts def = buildFromIdx(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);

std::mt19937 rng(sysSeed);
std::shuffle(masks.begin(), masks.end(), rng);
std::uniform_int_distribution<int> pick012(0, 2);

// Optional: keep unique combinations (by packing indices in 2 bits each)
auto packCode = [&](const std::array<int, 14>& idx) -> uint32_t {
uint32_t code = 0u;
for (int k = 0; k < 14; ++k) {
code |= (uint32_t(idx[k] & 0x3) << (2 * k));
}
return code;
};

sysCuts.clear();
sysCuts.reserve(1 + nSysRand);
sysCuts.push_back(def);
sysCuts.reserve(1 + (size_t)nSysRand);
sysCuts.push_back(def); // sysId=0

std::unordered_set<uint32_t> used;
used.reserve((size_t)nSysRand * 2);
used.insert(0u); // all-default code

const int nPick = std::max(0, (int)nSysRand);
while ((int)sysCuts.size() < 1 + nPick) {

std::array<int, 14> idx{};
for (int k = 0; k < 14; ++k)
idx[k] = pick012(rng);

uint32_t code = packCode(idx);
if (!used.insert(code).second)
continue; // already have this combination

const int nPick = std::min<int>(nSysRand, (int)masks.size());
for (int i = 0, picked = 0; picked < nPick && i < (int)masks.size(); ++i) {
if (masks[i] == 0u) { // avoid trivial mask that can reproduce default
// (optional) avoid generating exactly default again
if (code == 0u)
continue;
}
sysCuts.push_back(buildFromMask(masks[i]));
++picked;

sysCuts.push_back(buildFromIdx(
idx[0], idx[1], idx[2], idx[3],
idx[4], idx[5], idx[6], idx[7], idx[8], idx[9], idx[10],
idx[11], idx[12], idx[13]));
}

nSysTotal = (int)sysCuts.size();
}

Expand Down Expand Up @@ -598,7 +607,7 @@ struct f1protoncorrelation {
}
auto relative_momentum = getkstar(F1, Proton);
if (relative_momentum <= 0.5) {
histos.fill(HIST("hNsigmaProtonTPC"), protontrack.protonNsigmaTPC(), Proton.Pt());
histos.fill(HIST("hNsigmaProtonTPC"), protontrack.protonNsigmaTPC(), protontrack.protonNsigmaTOF(), Proton.Pt());
}
histos.fill(HIST("h2SameEventPtCorrelation"), relative_momentum, F1.Pt(), Proton.Pt());

Expand Down Expand Up @@ -968,7 +977,8 @@ struct f1protoncorrelation {
Kaon.SetXYZM(f1track.f1d2Px(), f1track.f1d2Py(), f1track.f1d2Pz(), 0.493);
Kshort.SetXYZM(f1track.f1d3Px(), f1track.f1d3Py(), f1track.f1d3Pz(), 0.497);
KaonKshortPair = Kaon + Kshort;

if (F1.Pt() < lowPtF1 || F1.Pt() > 50.0)
continue;
std::vector<int> activeSys;
activeSys.reserve((size_t)nSysTotal);

Expand Down Expand Up @@ -1014,7 +1024,7 @@ struct f1protoncorrelation {

const auto& sc0 = sysCuts[0];

if (countf1 && passPrimary(protontrack.protonDcaxy(), protontrack.protonDcaz(), protontrack.protonTPCNcrs(), protontrack.protonTPCNcls(), sc0)) {
if (countf1 && passPrimary(protontrack.protonDcaxy(), protontrack.protonDcaz(), protontrack.protonTPCNcrs(), protontrack.protonTPCNcls(), sc0) && passProtonPID(0, protontrack, Proton, pMinP, pMaxP, pTofP)) {
histos.fill(HIST("hNsigmaProtonTPC"), protontrack.protonNsigmaTPC(), protontrack.protonNsigmaTOF(), Proton.Pt());
}

Expand Down Expand Up @@ -1110,6 +1120,8 @@ struct f1protoncorrelation {
Kshort.SetXYZM(t1.f1d3Px(), t1.f1d3Py(), t1.f1d3Pz(), 0.497);
KaonKshortPair = Kaon + Kshort;
Proton.SetXYZM(t2.protonPx(), t2.protonPy(), t2.protonPz(), 0.938);
if (F1.Pt() < lowPtF1 || F1.Pt() > 50.0)
continue;
auto relative_momentum = getkstar(F1, Proton);
auto mT = getmT(F1, Proton);
// sys list for this (F1, p) pair
Expand Down
Loading
Loading