Skip to content

[PWGLF] Added effeciency calculation for polarisation study #18028

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
dmallick2 merged 1 commit into from
Sep 24, 2026
Merged

[PWGLF] Added effeciency calculation for polarisation study #18028

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
126 changes: 110 additions & 16 deletions PWGLF/Tasks/Resonances/lambda1520pbpb.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -269,21 +269,24 @@ struct Lambda1520pbpb {
histos.add("QAafter/Kaon/h2d_Kapr_nsigma_tof_p", " Kaons proton", kTH2F, {axisPpid, axisTOFNsigma});
histos.add("QAafter/Kaon/h2d_ka_nsigma_tof_vs_tpc", "n#sigma(TOF) vs n#sigma(TPC) Kaons", kTH2F, {axisTPCNsigma, axisTOFNsigma});

// Analysis
if (!doprocessMC) {
int nActiveFrames = static_cast<int>(cActivateHelicityFrame) +
static_cast<int>(cActivateCollinsSoperFrame) +
static_cast<int>(cActivateProductionFrame) +
static_cast<int>(cActivateBeamAxisFrame) +
static_cast<int>(cActivateRandomFrame);

bool polActive = nActiveFrames > 0;
// ── Polarization-frame activation (shared by data and MC booking) ───────
int nActiveFrames = static_cast<int>(cActivateHelicityFrame) +
static_cast<int>(cActivateCollinsSoperFrame) +
static_cast<int>(cActivateProductionFrame) +
static_cast<int>(cActivateBeamAxisFrame) +
static_cast<int>(cActivateRandomFrame);

bool polActive = nActiveFrames > 0;

if (nActiveFrames > 1) {
LOG(fatal) << "Multiple polarization frames activated, but axisPolFrame was removed! "
<< "Please activate only ONE frame at a time to prevent data mixing in the 4D sparse.";
}

if (nActiveFrames > 1) {
LOG(fatal) << "Multiple polarization frames activated, but axisPolFrame was removed! "
<< "Please activate only ONE frame at a time to prevent data mixing in the 4D sparse.";
}
const AxisSpec axisCosTheta{cCosThetaBins, "cos#theta*"};

// Analysis
if (!doprocessMC) {
if (!polActive) {
histos.add("Analysis/h4d_lstar_invm_US_PM", "THn #Lambda(1520)", kTHnSparseF, {axisInvM, axisPt, axisCent, axisOccupancy});
histos.add("Analysis/h4d_lstar_invm_US_MP", "THn #bar #Lambda(1520)", kTHnSparseF, {axisInvM, axisPt, axisCent, axisOccupancy});
Expand All @@ -305,8 +308,6 @@ struct Lambda1520pbpb {

// ── 4D Polarization / cosThetaStar histograms ──────────────────────────
if (polActive) {
const AxisSpec axisCosTheta{cCosThetaBins, "cos#theta*"};

histos.add("Analysis/h4d_pol_US_PM", "cos#theta* #Lambda(1520)", kTHnSparseF, {axisInvM, axisPt, axisCent, axisCosTheta});
histos.add("Analysis/h4d_pol_US_MP", "cos#theta* #bar{#Lambda}(1520)", kTHnSparseF, {axisInvM, axisPt, axisCent, axisCosTheta});
histos.add("Analysis/h4d_pol_LS_PP", "cos#theta* Like Signs p K^{+}", kTHnSparseF, {axisInvM, axisPt, axisCent, axisCosTheta});
Expand Down Expand Up @@ -340,6 +341,20 @@ struct Lambda1520pbpb {
histos.add("Analysis/h3d_rec_lstar_MP", "Reconstructed #bar{#Lambda}(1520) p_{T}", kTHnSparseF, {axisInvM, axisPt, axisCent});
histos.add("Analysis/h3d_reso_lstar_PM", "Resolution #Lambda(1520) p_{T}", kTHnSparseF, {{200, -0.05, 0.05}, axisPt, axisCent});
histos.add("Analysis/h3d_reso_lstar_MP", "Resolution #bar{#Lambda}(1520) p_{T}", kTHnSparseF, {{200, -0.05, 0.05}, axisPt, axisCent});

if (polActive) {
const AxisSpec axisMassRes{200, -0.05, 0.05, "M_{rec} - M_{gen} (GeV/#it{c}^{2})"};

histos.add("Analysis/h4d_pol_rec_PM", "cos#theta* Reconstructed (truth-matched) #Lambda(1520)", kTHnSparseF, {axisInvM, axisPt, axisCent, axisCosTheta});
histos.add("Analysis/h4d_pol_rec_MP", "cos#theta* Reconstructed (truth-matched) #bar{#Lambda}(1520)", kTHnSparseF, {axisInvM, axisPt, axisCent, axisCosTheta});
histos.add("Analysis/h4d_pol_reso_PM", "cos#theta* vs mass resolution #Lambda(1520)", kTHnSparseF, {axisMassRes, axisPt, axisCent, axisCosTheta});
histos.add("Analysis/h4d_pol_reso_MP", "cos#theta* vs mass resolution #bar{#Lambda}(1520)", kTHnSparseF, {axisMassRes, axisPt, axisCent, axisCosTheta});
}
}

if (doprocessMCGenPol && polActive) {
histos.add("Analysis/h4d_pol_gen_PM", "cos#theta* Generated (truth) #Lambda(1520)", kTHnSparseF, {axisInvM, axisPt, axisCent, axisCosTheta});
histos.add("Analysis/h4d_pol_gen_MP", "cos#theta* Generated (truth) #bar{#Lambda}(1520)", kTHnSparseF, {axisInvM, axisPt, axisCent, axisCosTheta});
}

if (doprocessMCGen) {
Expand Down Expand Up @@ -591,7 +606,10 @@ struct Lambda1520pbpb {
Rotated,
Mixed,
LikeSign,
LikeSignMixed };
LikeSignMixed,
McTrue,
McTrueRes,
GenTrue };

template <PolBkgMode Mode>
void fillPolarization(float candMass, float candPt,
Expand Down Expand Up @@ -629,6 +647,21 @@ struct Lambda1520pbpb {
histos.fill(HIST("Analysis/h4d_pol_LS_PP_mix"), candMass, candPt, mult, cosTheta);
else
histos.fill(HIST("Analysis/h4d_pol_LS_MM_mix"), candMass, candPt, mult, cosTheta);
} else if constexpr (Mode == PolBkgMode::McTrue) {
if (protonIsPositive)
histos.fill(HIST("Analysis/h4d_pol_rec_PM"), candMass, candPt, mult, cosTheta);
else
histos.fill(HIST("Analysis/h4d_pol_rec_MP"), candMass, candPt, mult, cosTheta);
} else if constexpr (Mode == PolBkgMode::McTrueRes) {
if (protonIsPositive)
histos.fill(HIST("Analysis/h4d_pol_reso_PM"), candMass, candPt, mult, cosTheta);
else
histos.fill(HIST("Analysis/h4d_pol_reso_MP"), candMass, candPt, mult, cosTheta);
} else if constexpr (Mode == PolBkgMode::GenTrue) {
if (protonIsPositive)
histos.fill(HIST("Analysis/h4d_pol_gen_PM"), candMass, candPt, mult, cosTheta);
else
histos.fill(HIST("Analysis/h4d_pol_gen_MP"), candMass, candPt, mult, cosTheta);
} else {
if (protonIsPositive)
histos.fill(HIST("Analysis/h4d_pol_US_PM"), candMass, candPt, mult, cosTheta);
Expand Down Expand Up @@ -975,6 +1008,17 @@ struct Lambda1520pbpb {
histos.fill(HIST("Analysis/h3d_rec_lstar_MP"), candMass, candPt, mult);
histos.fill(HIST("Analysis/h3d_reso_lstar_MP"), candMassRes, candPt, mult);
}

if (polActive) {
float eProton = std::sqrt(pxPr * pxPr + pyPr * pyPr + pzPr * pzPr + MassProton * MassProton);
float eKaon = std::sqrt(pxKa * pxKa + pyKa * pyKa + pzKa * pzKa + MassKaonCharged * MassKaonCharged);
std::array<float, 3> protonP = {pxPr, pyPr, pzPr};
std::array<float, 3> kaonP = {pxKa, pyKa, pzKa};
std::array<float, 3> motherP = {pxPr + pxKa, pyPr + pyKa, pzPr + pzKa};
float motherE = eProton + eKaon;
fillPolarization<PolBkgMode::McTrue>(candMass, candPt, motherP, motherE, protonP, eProton, kaonP, eKaon, mult, trkPr.motherPDG() > 0);
fillPolarization<PolBkgMode::McTrueRes>(candMassRes, candPt, motherP, motherE, protonP, eProton, kaonP, eKaon, mult, trkPr.motherPDG() > 0);
}
}
}
}
Expand Down Expand Up @@ -1073,6 +1117,56 @@ struct Lambda1520pbpb {
}
PROCESS_SWITCH(Lambda1520pbpb, processMC, "Process Event for MC", false);

void processMCGenPol(ResoMCCols::iterator const& collision, aod::ResoMCParents const& resoParents, aod::McParticles const& mcParticles)
{
auto mult = collision.cent();
bool polActive = static_cast<bool>(cActivateHelicityFrame) ||
static_cast<bool>(cActivateCollinsSoperFrame) ||
static_cast<bool>(cActivateProductionFrame) ||
static_cast<bool>(cActivateBeamAxisFrame) ||
static_cast<bool>(cActivateRandomFrame);
if (!polActive)
return;

for (auto const& part : resoParents) {
if (std::abs(part.pdgCode()) != lambda1520id)
continue;
float yshift = std::abs(part.y()) - cfgRapidityShift;
if (std::abs(yshift) > cfgRapidityCut)
continue;

bool pass1 = std::abs(part.daughterPDG1()) == kProton || std::abs(part.daughterPDG2()) == kProton;
bool pass2 = std::abs(part.daughterPDG1()) == kKPlus || std::abs(part.daughterPDG2()) == kKPlus;
if (!pass1 || !pass2)
continue;

std::array<float, 3> pvec = {part.px(), part.py(), part.pz()};
float mass = RecoDecay::m(pvec, part.e());

auto motherMc = mcParticles.iteratorAt(part.mcParticleId());
bool foundProton = false, foundKaon = false;
std::array<float, 3> protonP{}, kaonP{};
float protonE = 0.f, kaonE = 0.f;
for (auto const& dau : motherMc.daughters_as<aod::McParticles>()) {
if (!foundProton && std::abs(dau.pdgCode()) == kProton) {
protonP = {dau.px(), dau.py(), dau.pz()};
protonE = dau.e();
foundProton = true;
} else if (!foundKaon && std::abs(dau.pdgCode()) == kKPlus) {
kaonP = {dau.px(), dau.py(), dau.pz()};
kaonE = dau.e();
foundKaon = true;
}
}
if (!foundProton || !foundKaon)
continue;

std::array<float, 3> motherP = {part.px(), part.py(), part.pz()};
fillPolarization<PolBkgMode::GenTrue>(mass, part.pt(), motherP, part.e(), protonP, protonE, kaonP, kaonE, mult, part.pdgCode() > 0);
}
}
PROCESS_SWITCH(Lambda1520pbpb, processMCGenPol, "Standalone generator-level polarization pass (run separately from processMC)", false);

void processMCGen(ResoMCCols::iterator const& collision, aod::ResoMCParents const& resoParents)
{
float centrality = collision.cent();
Expand Down
Loading