Skip to content

Commit e5cfa50

Browse files
skundu692skundu692
authored
[PWGLF] Add new process function for MC spin correlation, fix systematic selection for f1-p correlation (#15230)
1 parent 8ff26f3 commit e5cfa50

File tree

4 files changed

+726
-89
lines changed

4 files changed

+726
-89
lines changed

‎PWGLF/TableProducer/Resonances/f1protonreducedtable.cxx‎

Lines changed: 52 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -84,6 +84,7 @@ struct f1protonreducedtable {
8484
Configurable<int> trackSphMin{"trackSphMin", 10, "Number of tracks for Spherocity Calculation"};
8585

8686
// Configs for track PID
87+
Configurable<bool> Confglobaltrackcheck{"Confglobaltrackcheck", true, "Global track check"};
8788
Configurable<bool> cfgSkimmedProcessing{"cfgSkimmedProcessing", true, "Analysed skimmed events"};
8889
Configurable<bool> ConfUseManualPIDproton{"ConfUseManualPIDproton", true, "True: use home-made PID solution for proton "};
8990
Configurable<bool> ConfUseManualPIDkaon{"ConfUseManualPIDkaon", true, "True: use home-made PID solution for kaon "};
@@ -161,16 +162,16 @@ struct f1protonreducedtable {
161162
{"hInvMassf1Like", "hInvMassf1Like", {HistType::kTH2F, {{400, 1.1f, 1.9f}, {100, 0.0f, 10.0f}}}},
162163
{"hInvMassf1kstar", "hInvMassf1kstar", {HistType::kTH3F, {{400, 1.1f, 1.9f}, {100, 0.0f, 10.0f}, {8, 0.0f, 0.8f}}}},
163164
{"hkstarDist", "hkstarDist", {HistType::kTH1F, {{300, 0.0f, 3.0f}}}},
164-
{"hDCAxy", "hDCAxy", {HistType::kTH1F, {{100, -5.0f, 5.0f}}}},
165-
{"hDCAz", "hDCAz", {HistType::kTH1F, {{100, -5.0f, 5.0f}}}},
165+
{"hDCAxy", "hDCAxy", {HistType::kTH3F, {{100, -0.05f, 0.05f}, {5, -2.5, 2.5}, {40, 0.0, 4.0}}}},
166+
{"hDCAz", "hDCAz", {HistType::kTH3F, {{100, -0.05f, 0.05f}, {2, 0, 2}, {40, 0.0, 4.0}}}},
166167
{"hPhi", "hPhi", {HistType::kTH1F, {{1400, -7.0f, 7.0f}}}},
167168
{"hPhiSphero", "hPhiSphero", {HistType::kTH1F, {{1400, -7.0f, 7.0f}}}},
168169
{"hEta", "hEta", {HistType::kTH1F, {{20, -1.0f, 1.0f}}}},
169170
{"hNsigmaPtpionTPC", "hNsigmaPtpionTPC", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
170171
{"hNsigmaPtpionTOF", "hNsigmaPtpionTOF", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
171172
{"hNsigmaPtkaonTPC", "hNsigmaPtkaonTPC", {HistType::kTH3F, {{200, -10.0f, 10.0f}, {200, -20.0f, 20.0f}, {100, 0.0f, 10.0f}}}},
172173
{"hNsigmaPtkaonTOF", "hNsigmaPtkaonTOF", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
173-
{"hNsigmaPtprotonTPC", "hNsigmaPtprotonTPC", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
174+
{"hNsigmaPtprotonTPC", "hNsigmaPtprotonTPC", {HistType::kTH3F, {{100, -5.0f, 5.0f}, {100, -5.0f, 5.0f}, {100, 0.0f, 10.0f}}}},
174175
{"hNsigmaPtprotonTOF", "hNsigmaPtprotonTOF", {HistType::kTH2F, {{200, -10.0f, 10.0f}, {100, 0.0f, 10.0f}}}},
175176
{"hInvMassk0", "hInvMassk0", {HistType::kTH2F, {{200, 0.4f, 0.6f}, {100, 0.0f, 10.0f}}}},
176177
},
@@ -265,7 +266,7 @@ struct f1protonreducedtable {
265266
template <typename T>
266267
bool selectionGlobalTrack(const T& candidate)
267268
{
268-
if (!(candidate.isGlobalTrack() && candidate.isPVContributor())) {
269+
if (Confglobaltrackcheck && !(candidate.isGlobalTrack() && candidate.isPVContributor())) {
269270
return false;
270271
}
271272
return true;
@@ -348,6 +349,32 @@ struct f1protonreducedtable {
348349
return false;
349350
}
350351

352+
inline bool passProtonPID(float nsigmaTPC, float nsigmaTOF, float TOFHit, ROOT::Math::PtEtaPhiMVector proton)
353+
{
354+
// pidMode:
355+
// 0 = default: p < thr -> |TPC| < 2.5 ; p >= thr -> TOF mandatory AND circular(TPC,TOF) < 2.0
356+
// 1 = syst-1: p < thr -> |TPC| < 2.0 ; p >= thr -> TOF mandatory AND circular(TPC,TOF) < 2.0
357+
// 2 = syst-2: p < thr -> |TPC| < 2.5 ; p >= thr -> TOF mandatory AND circular(TPC,TOF) < 2.5
358+
359+
if (proton.Pt() > 4.0 || proton.Pt() < 0.15) {
360+
return false;
361+
}
362+
363+
if (proton.P() < 0.7) {
364+
return (std::abs(nsigmaTPC) < 2.5);
365+
}
366+
367+
// above threshold: TOF must exist
368+
if (TOFHit != 1) {
369+
return false;
370+
}
371+
372+
const float nsTPC = nsigmaTPC;
373+
const float nsTOF = nsigmaTOF;
374+
const float comb = std::sqrt(nsTPC * nsTPC + nsTOF * nsTOF);
375+
return (comb < 2.5);
376+
}
377+
351378
template <typename Collision, typename V0>
352379
bool SelectionV0(Collision const& collision, V0 const& candidate)
353380
{
@@ -682,8 +709,6 @@ struct f1protonreducedtable {
682709
if (!selectionGlobalTrack(track)) {
683710
continue;
684711
}
685-
qaRegistry.fill(HIST("hDCAxy"), track.dcaXY());
686-
qaRegistry.fill(HIST("hDCAz"), track.dcaZ());
687712
qaRegistry.fill(HIST("hEta"), track.eta());
688713
qaRegistry.fill(HIST("hPhi"), track.phi());
689714
double nTPCSigmaP[3]{track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()};
@@ -778,17 +803,15 @@ struct f1protonreducedtable {
778803
ProtonIndex.push_back(track.globalIndex());
779804
ProtonCharge.push_back(track.sign());
780805

781-
ProtonDcaxy.push_back(std::abs(track.dcaXY()));
782-
ProtonDcaz.push_back(std::abs(track.dcaZ()));
806+
ProtonDcaxy.push_back(track.dcaXY());
807+
ProtonDcaz.push_back(track.dcaZ());
783808
ProtonTPCNcls.push_back(std::abs(track.tpcNClsFound()));
784809
ProtonTPCNcrs.push_back(std::abs(track.tpcNClsCrossedRows()));
785810

786811
if (track.sign() > 0) {
787-
qaRegistry.fill(HIST("hNsigmaPtprotonTPC"), nTPCSigmaP[2], track.pt());
788812
ProtonTPCNsigma.push_back(nTPCSigmaP[2]);
789813
}
790814
if (track.sign() < 0) {
791-
qaRegistry.fill(HIST("hNsigmaPtprotonTPC"), nTPCSigmaN[2], track.pt());
792815
ProtonTPCNsigma.push_back(nTPCSigmaN[2]);
793816
}
794817
if (track.hasTOF()) {
@@ -961,6 +984,25 @@ struct f1protonreducedtable {
961984
}
962985
}
963986
qaRegistry.fill(HIST("hEventstat"), 0.5);
987+
if (keepEventF1Proton) {
988+
for (auto iproton = protons.begin(); iproton != protons.end(); ++iproton) {
989+
auto i6 = std::distance(protons.begin(), iproton);
990+
ProtonVectorDummy2 = protons.at(i6);
991+
if (std::abs(ProtonDcaxy.at(i6)) < 0.05 && std::abs(ProtonDcaz.at(i6)) < 0.05) {
992+
if (ProtonTOFHit.at(i6) && ProtonVectorDummy2.P() > 0.7) {
993+
qaRegistry.fill(HIST("hNsigmaPtprotonTPC"), ProtonTPCNsigma.at(i6), ProtonTOFNsigma.at(i6), ProtonVectorDummy2.Pt());
994+
}
995+
if (ProtonVectorDummy2.P() < 0.7) {
996+
qaRegistry.fill(HIST("hNsigmaPtprotonTPC"), ProtonTPCNsigma.at(i6), 4.999, ProtonVectorDummy2.Pt());
997+
}
998+
}
999+
if (passProtonPID(ProtonTPCNsigma.at(i6), ProtonTOFNsigma.at(i6), ProtonTOFHit.at(i6), ProtonVectorDummy2)) {
1000+
qaRegistry.fill(HIST("hDCAxy"), ProtonDcaxy.at(i6), ProtonCharge.at(i6), ProtonVectorDummy2.Pt());
1001+
qaRegistry.fill(HIST("hDCAz"), ProtonDcaz.at(i6), ProtonCharge.at(i6), ProtonVectorDummy2.Pt());
1002+
}
1003+
}
1004+
}
1005+
9641006
if (numberF1 > 0 && (f1resonance.size() == f1signal.size()) && (f1resonance.size() == f1kaonkshortmass.size()) && (f1resonance.size() == f1resonanced1.size()) && (f1resonance.size() == f1resonanced2.size()) && (f1resonance.size() == f1resonanced3.size())) {
9651007
qaRegistry.fill(HIST("hEventstat"), 1.5);
9661008
if (keepEventF1Proton) {

‎PWGLF/Tasks/Resonances/f1protoncorrelation.cxx‎

Lines changed: 84 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,7 @@
4444
#include <random>
4545
#include <sstream>
4646
#include <string>
47+
#include <unordered_set>
4748
#include <vector>
4849

4950
using namespace o2;
@@ -140,86 +141,94 @@ struct f1protoncorrelation {
140141
// -------------------------
141142
void buildSystematicCuts()
142143
{
143-
// choices from your table/picture
144-
const std::array<float, 2> optDcaxy{0.01f, 0.03f};
145-
const std::array<float, 2> optDcaz{0.01f, 0.03f};
146-
const std::array<int, 2> optNcrs{90, 100};
147-
const std::array<int, 2> optNcls{90, 100};
148-
149-
const std::array<float, 2> optCPA{0.99f, 0.995f};
150-
const std::array<float, 2> optRad{0.8f, 1.0f};
151-
const std::array<float, 2> optDcaDD{0.9f, 0.8f};
152-
const std::array<float, 2> optDcaV0{0.2f, 0.15f};
153-
const std::array<float, 2> optLife{16.f, 18.f};
154-
const std::array<float, 2> optDcaD1{0.06f, 0.08f};
155-
const std::array<float, 2> optDcaD2{0.06f, 0.08f};
156-
157-
// bit layout:
158-
// primary: 4 bits (0..3)
159-
// v0: 7 bits (4..10)
160-
// pid: 3 bits (11..13)
161-
// total: 14 bits -> 16384 combos
162-
auto buildFromMask = [&](uint32_t m) -> SysCuts {
144+
// 3 options per cut: index 0 = DEFAULT, index 1/2 = variations
145+
// Fill these with the exact values you want (I used your def as index 0 + your old options as 1/2)
146+
const std::array<float, 3> optDcaxy{0.05f, 0.01f, 0.03f};
147+
const std::array<float, 3> optDcaz{0.05f, 0.01f, 0.03f};
148+
const std::array<int, 3> optNcrs{80, 90, 100};
149+
const std::array<int, 3> optNcls{80, 90, 100};
150+
151+
const std::array<float, 3> optCPA{0.985f, 0.99f, 0.995f};
152+
const std::array<float, 3> optRad{0.50f, 0.80f, 1.00f};
153+
const std::array<float, 3> optDcaDD{1.00f, 0.90f, 0.80f};
154+
const std::array<float, 3> optDcaV0{0.30f, 0.20f, 0.15f};
155+
const std::array<float, 3> optLife{20.f, 16.f, 18.f};
156+
const std::array<float, 3> optDcaD1{0.05f, 0.06f, 0.08f};
157+
const std::array<float, 3> optDcaD2{0.05f, 0.06f, 0.08f};
158+
159+
// PID modes: also allow default (0) to appear in random variations
160+
const std::array<int, 3> optPidPi{0, 1, 2};
161+
const std::array<int, 3> optPidK{0, 1, 2};
162+
const std::array<int, 3> optPidP{0, 1, 2};
163+
164+
// Helper: build SysCuts from chosen indices (0..2)
165+
auto buildFromIdx = [&](int i0, int i1, int i2, int i3,
166+
int i4, int i5, int i6, int i7, int i8, int i9, int i10,
167+
int i11, int i12, int i13) -> SysCuts {
163168
SysCuts c{};
164-
c.maxDcaxy = optDcaxy[(m >> 0) & 1u];
165-
c.maxDcaz = optDcaz[(m >> 1) & 1u];
166-
c.minTPCCrossedRows = optNcrs[(m >> 2) & 1u];
167-
c.minTPCClusters = optNcls[(m >> 3) & 1u];
168-
169-
c.minCPA = optCPA[(m >> 4) & 1u];
170-
c.minRadius = optRad[(m >> 5) & 1u];
171-
c.maxDcaDaughters = optDcaDD[(m >> 6) & 1u];
172-
c.maxDcaV0 = optDcaV0[(m >> 7) & 1u];
173-
c.maxLifetime = optLife[(m >> 8) & 1u];
174-
c.minDcaD1 = optDcaD1[(m >> 9) & 1u];
175-
c.minDcaD2 = optDcaD2[(m >> 10) & 1u];
176-
177-
c.pidPi = (m >> 11) & 1u;
178-
c.pidK = (m >> 12) & 1u;
179-
c.pidP = (m >> 13) & 1u;
169+
c.maxDcaxy = optDcaxy[i0];
170+
c.maxDcaz = optDcaz[i1];
171+
c.minTPCCrossedRows = optNcrs[i2];
172+
c.minTPCClusters = optNcls[i3];
173+
174+
c.minCPA = optCPA[i4];
175+
c.minRadius = optRad[i5];
176+
c.maxDcaDaughters = optDcaDD[i6];
177+
c.maxDcaV0 = optDcaV0[i7];
178+
c.maxLifetime = optLife[i8];
179+
c.minDcaD1 = optDcaD1[i9];
180+
c.minDcaD2 = optDcaD2[i10];
181+
182+
c.pidPi = optPidPi[i11];
183+
c.pidK = optPidK[i12];
184+
c.pidP = optPidP[i13];
180185
return c;
181186
};
182187

183-
// DEFAULT (sysId=0): set to your baseline
184-
SysCuts def{};
185-
def.maxDcaxy = 0.05f;
186-
def.maxDcaz = 0.05f;
187-
def.minTPCCrossedRows = 80;
188-
def.minTPCClusters = 80;
189-
190-
def.minCPA = 0.985f;
191-
def.minRadius = 0.5f;
192-
def.maxDcaDaughters = 1.0f;
193-
def.maxDcaV0 = 0.3f;
194-
def.maxLifetime = 20.f;
195-
def.minDcaD1 = 0.05f;
196-
def.minDcaD2 = 0.05f;
197-
198-
def.pidPi = 0;
199-
def.pidK = 0;
200-
def.pidP = 0;
201-
202-
std::vector<uint32_t> masks;
203-
masks.reserve(16384);
204-
for (uint32_t m = 0; m < 16384; ++m) {
205-
masks.push_back(m);
206-
}
188+
// sysId=0 must be strict default (all indices = 0)
189+
SysCuts def = buildFromIdx(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
207190

208191
std::mt19937 rng(sysSeed);
209-
std::shuffle(masks.begin(), masks.end(), rng);
192+
std::uniform_int_distribution<int> pick012(0, 2);
193+
194+
// Optional: keep unique combinations (by packing indices in 2 bits each)
195+
auto packCode = [&](const std::array<int, 14>& idx) -> uint32_t {
196+
uint32_t code = 0u;
197+
for (int k = 0; k < 14; ++k) {
198+
code |= (uint32_t(idx[k] & 0x3) << (2 * k));
199+
}
200+
return code;
201+
};
210202

211203
sysCuts.clear();
212-
sysCuts.reserve(1 + nSysRand);
213-
sysCuts.push_back(def);
204+
sysCuts.reserve(1 + (size_t)nSysRand);
205+
sysCuts.push_back(def); // sysId=0
206+
207+
std::unordered_set<uint32_t> used;
208+
used.reserve((size_t)nSysRand * 2);
209+
used.insert(0u); // all-default code
210+
211+
const int nPick = std::max(0, (int)nSysRand);
212+
while ((int)sysCuts.size() < 1 + nPick) {
213+
214+
std::array<int, 14> idx{};
215+
for (int k = 0; k < 14; ++k)
216+
idx[k] = pick012(rng);
217+
218+
uint32_t code = packCode(idx);
219+
if (!used.insert(code).second)
220+
continue; // already have this combination
214221

215-
const int nPick = std::min<int>(nSysRand, (int)masks.size());
216-
for (int i = 0, picked = 0; picked < nPick && i < (int)masks.size(); ++i) {
217-
if (masks[i] == 0u) { // avoid trivial mask that can reproduce default
222+
// (optional) avoid generating exactly default again
223+
if (code == 0u)
218224
continue;
219-
}
220-
sysCuts.push_back(buildFromMask(masks[i]));
221-
++picked;
225+
226+
sysCuts.push_back(buildFromIdx(
227+
idx[0], idx[1], idx[2], idx[3],
228+
idx[4], idx[5], idx[6], idx[7], idx[8], idx[9], idx[10],
229+
idx[11], idx[12], idx[13]));
222230
}
231+
223232
nSysTotal = (int)sysCuts.size();
224233
}
225234

@@ -598,7 +607,7 @@ struct f1protoncorrelation {
598607
}
599608
auto relative_momentum = getkstar(F1, Proton);
600609
if (relative_momentum <= 0.5) {
601-
histos.fill(HIST("hNsigmaProtonTPC"), protontrack.protonNsigmaTPC(), Proton.Pt());
610+
histos.fill(HIST("hNsigmaProtonTPC"), protontrack.protonNsigmaTPC(), protontrack.protonNsigmaTOF(), Proton.Pt());
602611
}
603612
histos.fill(HIST("h2SameEventPtCorrelation"), relative_momentum, F1.Pt(), Proton.Pt());
604613

@@ -968,7 +977,8 @@ struct f1protoncorrelation {
968977
Kaon.SetXYZM(f1track.f1d2Px(), f1track.f1d2Py(), f1track.f1d2Pz(), 0.493);
969978
Kshort.SetXYZM(f1track.f1d3Px(), f1track.f1d3Py(), f1track.f1d3Pz(), 0.497);
970979
KaonKshortPair = Kaon + Kshort;
971-
980+
if (F1.Pt() < lowPtF1 || F1.Pt() > 50.0)
981+
continue;
972982
std::vector<int> activeSys;
973983
activeSys.reserve((size_t)nSysTotal);
974984

@@ -1014,7 +1024,7 @@ struct f1protoncorrelation {
10141024

10151025
const auto& sc0 = sysCuts[0];
10161026

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

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

0 commit comments

Comments
 (0)