Skip to content

[PWGLF] added DCAxy p-dependent cut and added rotational method #15274

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
ktf merged 2 commits into from
Mar 6, 2026
+127 −37
Merged

[PWGLF] added DCAxy p-dependent cut and added rotational method #15274

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3d952f3
added p dep. dcaxy cut and updated rota bg method
nasirmehdimalik Mar 5, 2026
7d54e04
Merge branch 'AliceO2Group:master' into master
nasirmehdimalik Mar 5, 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
164 changes: 127 additions & 37 deletions PWGLF/Tasks/Resonances/lambda1520_PbPb.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,9 +29,6 @@
#include <Framework/HistogramRegistry.h>
#include <Framework/HistogramSpec.h>

#include <TLorentzVector.h>
#include <TRandom.h>

#include <fairlogger/Logger.h>

using namespace o2;
Expand All @@ -56,7 +53,15 @@ struct lambdaAnalysis_pb {
Configurable<float> cPMin{"cPMin", 0., "Minimum Track p"};
Configurable<float> cEtaCut{"cEtaCut", 0.8, "Pseudorapidity cut"};
Configurable<float> cDcaz{"cDcazMin", 1., "Minimum DCAz"};
Configurable<float> cDcaxy{"cDcaxyMin", 0.1, "Minimum DCAxy"};

Configurable<std::vector<float>> cDcaPtBinsPr{"cDcaPtBinsPr", {0.0f, 0.5f, 1.0f, 2.0f, 3.0f, 5.0f, 1000.0f}, "Proton pT bin edges for DCAxy cut"};

Configurable<std::vector<float>> cDcaXYBinsPr{"cDcaXYBinsPr", {0.020f, 0.015f, 0.010f, 0.007f, 0.005f, 0.004f}, "Proton max |DCAxy| per pT bin (cm)"};

// Kaon DCAxy — pT binned
Configurable<std::vector<float>> cDcaPtBinsKa{"cDcaPtBinsKa", {0.0f, 0.3f, 0.6f, 1.0f, 2.0f, 1000.0f}, "Kaon pT bin edges for DCAxy cut"};

Configurable<std::vector<float>> cDcaXYBinsKa{"cDcaXYBinsKa", {0.025f, 0.018f, 0.012f, 0.008f, 0.004f}, "Kaon max |DCAxy| per pT bin (cm)"};
Configurable<bool> isonlyQC{"isonlyQC", false, "only QC"};
Configurable<bool> isDeepAngle{"isDeepAngle", false, "Deep Angle cut"};
Configurable<double> cfgDeepAngle{"cfgDeepAngle", 0.04, "Deep Angle cut value"};
Expand Down Expand Up @@ -108,6 +113,16 @@ struct lambdaAnalysis_pb {
ConfigurableAxis cMixMultBins{"cMixMultBins", {VARIABLE_WIDTH, 0.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 200.0f}, "Mixing bins - multiplicity"};
ConfigurableAxis cMixEPAngle{"cMixEPAngle", {VARIABLE_WIDTH, -1.5708f, -1.25664f, -0.942478f, -0.628319f, 0.f, 0.628319f, 0.942478f, 1.25664f, 1.5708f}, "event plane"};
ConfigurableAxis occupancy_bins{"occupancy_bins", {VARIABLE_WIDTH, 0.0, 100, 500, 600, 1000, 1100, 1500, 1600, 2000, 2100, 2500, 2600, 3000, 3100, 3500, 3600, 4000, 4100, 4500, 4600, 5000, 5100, 9999}, "Binning of the occupancy axis"};
Configurable<int> cNofRotations{"cNofRotations", 10, "Number of rotations for rotational background"};
Configurable<float> rotationalcut{"rotationalcut", 6.f, "Rotational background angle window: PI/rotationalcut"};

// ── MC Event Selection Configurables ─────────────────────────────────────
Configurable<bool> cEvtMCAfterAllCuts{"cEvtMCAfterAllCuts", false, "MC event sel: isInAfterAllCuts"};
Configurable<bool> cEvtMCINELgt0{"cEvtMCINELgt0", false, "MC event sel: isINELgt0"};
Configurable<bool> cEvtMCSel8{"cEvtMCSel8", false, "MC event sel: isInSel8"};
Configurable<bool> cEvtMCVtxIn10{"cEvtMCVtxIn10", false, "MC event sel: isVtxIn10"};
Configurable<bool> cEvtMCTriggerTVX{"cEvtMCTriggerTVX", false, "MC event sel: isTriggerTVX"};
Configurable<bool> cEvtMCRecINELgt0{"cEvtMCRecINELgt0", false, "MC event sel: isRecINELgt0"};
// Histogram Registry.
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};

Expand All @@ -131,6 +146,7 @@ struct lambdaAnalysis_pb {
AxisSpec axisDCAz = {cDCAzBins, "DCA_{z} (cm)"};

histos.add("Event/h1d_ft0_mult_percentile", "FT0 (%)", kTH2F, {axisCent, axisOccupancy});
histos.add("Event/h_ft0_vz", "Collision Vertex Z position", kTH1F, {{100, -15., 15.}});
if (doprocessMix || doprocessMixDF || doprocessMixepDF) {
histos.add("Event/mixing_vzVsmultpercentile", "FT0(%)", kTH3F, {axisCent, axisVz, axisEP});
}
Expand Down Expand Up @@ -208,19 +224,58 @@ struct lambdaAnalysis_pb {
if (std::abs(track.eta()) > cEtaCut)
return false;

if (std::abs(track.dcaZ()) > cDcaz)
if (cPrimaryTrack && !track.isPrimaryTrack())
return false;

if (std::abs(track.dcaXY()) > cDcaxy)
if (cGlobalWoDCATrack && !track.isGlobalTrackWoDCA())
return false;

if (cPrimaryTrack && !track.isPrimaryTrack())
if (cPVContributor && !track.isPVContributor())
return false;

if (cGlobalWoDCATrack && !track.isGlobalTrackWoDCA())
return true;
}
// ── Proton DCA Selection ──────────────────────────────────────────────────
template <typename T>
bool dcaSelectionProton(T const& track, float p)
{
auto ptBinsPr = static_cast<std::vector<float>>(cDcaPtBinsPr);
auto dcaXYPr = static_cast<std::vector<float>>(cDcaXYBinsPr);
int nBinsPr = static_cast<int>(ptBinsPr.size()) - 1;

bool dcaXYPassed = false;
for (int i = 0; i < nBinsPr; i++) {
if (p >= ptBinsPr[i] && p < ptBinsPr[i + 1] &&
std::abs(track.dcaXY()) < dcaXYPr[i])
dcaXYPassed = true;
}
if (!dcaXYPassed)
return false;

if (cPVContributor && !track.isPVContributor())
if (std::abs(track.dcaZ()) > cDcaz)
return false;

return true;
}

// ── Kaon DCA Selection ────────────────────────────────────────────────────
template <typename T>
bool dcaSelectionKaon(T const& track, float p)
{
auto ptBinsKa = static_cast<std::vector<float>>(cDcaPtBinsKa);
auto dcaXYKa = static_cast<std::vector<float>>(cDcaXYBinsKa);
int nBinsKa = static_cast<int>(ptBinsKa.size()) - 1;

bool dcaXYPassed = false;
for (int i = 0; i < nBinsKa; i++) {
if (p >= ptBinsKa[i] && p < ptBinsKa[i + 1] &&
std::abs(track.dcaXY()) < dcaXYKa[i])
dcaXYPassed = true;
}
if (!dcaXYPassed)
return false;

if (std::abs(track.dcaZ()) > cDcaz)
return false;

return true;
Expand Down Expand Up @@ -362,8 +417,6 @@ struct lambdaAnalysis_pb {
void fillDataHistos(trackType const& trk1, trackType const& trk2, float mult, int occup = 100)
{

TLorentzVector p1, p2, p;
TRandom* rn = new TRandom();
float p_ptot = 0., k_ptot = 0.;

for (auto const& [trkPr, trkKa] : soa::combinations(soa::CombinationsFullIndexPolicy(trk1, trk2))) {
Expand Down Expand Up @@ -412,6 +465,9 @@ struct lambdaAnalysis_pb {
continue;
if (!selectionPIDProton(trkPr, p_ptot) || !selectionPIDKaon(trkKa, k_ptot))
continue;
if (!dcaSelectionProton(trkPr, p_ptot) || !dcaSelectionKaon(trkKa, k_ptot))
continue;

if (isDeepAngle && TMath::ACos((trkPr.pt() * trkKa.pt() + _pzPr * _pzKa) / (p_ptot * k_ptot)) < cfgDeepAngle)
continue;

Expand Down Expand Up @@ -461,39 +517,56 @@ struct lambdaAnalysis_pb {

if (isonlyQC)
continue;
// Invariant mass reconstruction.
p1.SetXYZM(_pxPr, _pyPr, _pzPr, MassProton);
p2.SetXYZM(_pxKa, _pyKa, _pzKa, MassKaonCharged);
p = p1 + p2;

if (std::abs(p.Rapidity()) > 0.5)
continue;
std::array<float, 3> pvec0 = {_pxPr, _pyPr, _pzPr};
std::array<float, 3> pvec1 = {_pxKa, _pyKa, _pzKa};
std::array<std::array<float, 3>, 2> arrMomrec = {pvec0, pvec1};

float _M = RecoDecay::m(arrMomrec, std::array{MassProton, MassKaonCharged});
float _pt = RecoDecay::pt(std::array{_pxPr + _pxKa, _pyPr + _pyKa});

auto _M = p.M();
auto _pt = p.Pt();
if (std::abs(RecoDecay::y(std::array{_pxPr + _pxKa, _pyPr + _pyKa, _pzPr + _pzKa}, MassLambda1520)) > 0.5)
continue;

// Apply kinematic cuts.
if (cKinCuts) {
TVector3 v1(_pxPr, _pyPr, _pzPr);
TVector3 v2(_pxKa, _pyKa, _pzKa);
float alpha = v1.Angle(v2);
if (alpha > 1.4 && alpha < 2.4)
continue;
}

// Fill Invariant Mass Histograms.
if constexpr (!mix && !mc) {
if (trkPr.sign() * trkKa.sign() < 0) {

if (trkPr.sign() > 0)
histos.fill(HIST("Analysis/h4d_lstar_invm_US_PM"), _M, _pt, mult, occup);
else
histos.fill(HIST("Analysis/h4d_lstar_invm_US_MP"), _M, _pt, mult, occup);

if (doRotate) {
float theta = rn->Uniform(1.56, 1.58);
p1.RotateZ(theta);
p = p1 + p2;
if (std::abs(p.Rapidity()) < 0.5) {
histos.fill(HIST("Analysis/h4d_lstar_invm_rot"), p.M(), p.Pt(), mult, occup);
for (int i = 0; i < cNofRotations; i++) {

float delta = o2::constants::math::PI / rotationalcut;
float theta2 = (o2::constants::math::PI - delta) + i * (2.f * delta / (cNofRotations - 1));

float phiRot = trkKa.phi() + theta2;
if (phiRot > o2::constants::math::TwoPI)
phiRot -= o2::constants::math::TwoPI;
if (phiRot < 0.f)
phiRot += o2::constants::math::TwoPI;

float _pxKaRot = trkKa.pt() * std::cos(phiRot);
float _pyKaRot = trkKa.pt() * std::sin(phiRot);

std::array<float, 3> pvec0rot = {_pxPr, _pyPr, _pzPr};
std::array<float, 3> pvec1rot = {_pxKaRot, _pyKaRot, _pzKa};
std::array<std::array<float, 3>, 2> arrMomRot = {pvec0rot, pvec1rot};

float _Mrot = RecoDecay::m(arrMomRot, std::array{MassProton, MassKaonCharged});
float _ptRot = RecoDecay::pt(std::array{_pxPr + _pxKaRot, _pyPr + _pyKaRot});

if (std::abs(RecoDecay::y(
std::array{_pxPr + _pxKaRot, _pyPr + _pyKaRot, _pzPr + _pzKa},
MassLambda1520)) > 0.5f)
continue;

histos.fill(HIST("Analysis/h4d_lstar_invm_rot"), _Mrot, _ptRot, mult, occup);
}
}
} else {
Expand Down Expand Up @@ -542,24 +615,41 @@ struct lambdaAnalysis_pb {
}

using resoCols = soa::Join<aod::ResoCollisions, aod::ResoEvtPlCollisions>;
using resoMCCols = soa::Join<aod::ResoCollisions, aod::ResoMCCollisions>;

using resoTracks = aod::ResoTracks;

void processData(resoCols::iterator const& collision, resoTracks const& tracks)
{

// LOGF(info, " collisions: Index = %d %d", collision.globalIndex(),tracks.size());
histos.fill(HIST("Event/h1d_ft0_mult_percentile"), collision.cent(), 100);
histos.fill(HIST("Event/h_ft0_vz"), collision.posZ());

fillDataHistos<false, false>(tracks, tracks, collision.cent());
}

PROCESS_SWITCH(lambdaAnalysis_pb, processData, "Process for Same Event Data", true);

void processMC(resoCols::iterator const& collision,
void processMC(resoMCCols::iterator const& collision,
soa::Join<aod::ResoTracks, aod::ResoMCTracks> const& tracks, aod::ResoMCParents const& resoParents)
{
if (cEvtMCAfterAllCuts && !collision.isInAfterAllCuts())
return;
if (cEvtMCINELgt0 && !collision.isINELgt0())
return;
if (cEvtMCSel8 && !collision.isInSel8())
return;
if (cEvtMCVtxIn10 && !collision.isVtxIn10())
return;
if (cEvtMCTriggerTVX && !collision.isTriggerTVX())
return;
if (cEvtMCRecINELgt0 && !collision.isRecINELgt0())
return;

auto mult = collision.cent();
histos.fill(HIST("Event/h1d_ft0_mult_percentile"), mult, 100);
histos.fill(HIST("Event/h_ft0_vz"), collision.posZ());
fillDataHistos<false, true>(tracks, tracks, mult);

// get MC pT-spectra
Expand Down Expand Up @@ -587,9 +677,9 @@ struct lambdaAnalysis_pb {
}

for (auto const& part : resoParents) {
if (abs(part.pdgCode()) != lambda1520id) // // L* pdg_code = 3124
if (std::abs(part.pdgCode()) != lambda1520id) // // L* pdg_code = 3124
continue;
if (abs(part.y()) > 0.5) { // rapidity cut
if (std::abs(part.y()) > 0.5) { // rapidity cut
continue;
}

Expand All @@ -605,10 +695,10 @@ struct lambdaAnalysis_pb {

if (!pass1 || !pass2) // If we have both decay products
continue;
std::array<float, 3> pvec = {part.px(), part.py(), part.pz()};

float mass = RecoDecay::m(pvec, part.e());

TLorentzVector p4;
p4.SetPxPyPzE(part.px(), part.py(), part.pz(), part.e());
auto mass = p4.M();
if (part.pdgCode() > 0)
histos.fill(HIST("Analysis/h3d_gen_lstar_PM"), mass, part.pt(), mult);
else
Expand Down
Loading