diff --git a/PWGLF/Tasks/Resonances/lambda1520_PbPb.cxx b/PWGLF/Tasks/Resonances/lambda1520_PbPb.cxx index d8edd64cdb9..6af3bde09f8 100644 --- a/PWGLF/Tasks/Resonances/lambda1520_PbPb.cxx +++ b/PWGLF/Tasks/Resonances/lambda1520_PbPb.cxx @@ -29,9 +29,6 @@ #include #include -#include -#include - #include using namespace o2; @@ -56,7 +53,15 @@ struct lambdaAnalysis_pb { Configurable cPMin{"cPMin", 0., "Minimum Track p"}; Configurable cEtaCut{"cEtaCut", 0.8, "Pseudorapidity cut"}; Configurable cDcaz{"cDcazMin", 1., "Minimum DCAz"}; - Configurable cDcaxy{"cDcaxyMin", 0.1, "Minimum DCAxy"}; + + Configurable> cDcaPtBinsPr{"cDcaPtBinsPr", {0.0f, 0.5f, 1.0f, 2.0f, 3.0f, 5.0f, 1000.0f}, "Proton pT bin edges for DCAxy cut"}; + + Configurable> 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> cDcaPtBinsKa{"cDcaPtBinsKa", {0.0f, 0.3f, 0.6f, 1.0f, 2.0f, 1000.0f}, "Kaon pT bin edges for DCAxy cut"}; + + Configurable> cDcaXYBinsKa{"cDcaXYBinsKa", {0.025f, 0.018f, 0.012f, 0.008f, 0.004f}, "Kaon max |DCAxy| per pT bin (cm)"}; Configurable isonlyQC{"isonlyQC", false, "only QC"}; Configurable isDeepAngle{"isDeepAngle", false, "Deep Angle cut"}; Configurable cfgDeepAngle{"cfgDeepAngle", 0.04, "Deep Angle cut value"}; @@ -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 cNofRotations{"cNofRotations", 10, "Number of rotations for rotational background"}; + Configurable rotationalcut{"rotationalcut", 6.f, "Rotational background angle window: PI/rotationalcut"}; + + // ── MC Event Selection Configurables ───────────────────────────────────── + Configurable cEvtMCAfterAllCuts{"cEvtMCAfterAllCuts", false, "MC event sel: isInAfterAllCuts"}; + Configurable cEvtMCINELgt0{"cEvtMCINELgt0", false, "MC event sel: isINELgt0"}; + Configurable cEvtMCSel8{"cEvtMCSel8", false, "MC event sel: isInSel8"}; + Configurable cEvtMCVtxIn10{"cEvtMCVtxIn10", false, "MC event sel: isVtxIn10"}; + Configurable cEvtMCTriggerTVX{"cEvtMCTriggerTVX", false, "MC event sel: isTriggerTVX"}; + Configurable cEvtMCRecINELgt0{"cEvtMCRecINELgt0", false, "MC event sel: isRecINELgt0"}; // Histogram Registry. HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -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}); } @@ -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 + bool dcaSelectionProton(T const& track, float p) + { + auto ptBinsPr = static_cast>(cDcaPtBinsPr); + auto dcaXYPr = static_cast>(cDcaXYBinsPr); + int nBinsPr = static_cast(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 + bool dcaSelectionKaon(T const& track, float p) + { + auto ptBinsKa = static_cast>(cDcaPtBinsKa); + auto dcaXYKa = static_cast>(cDcaXYBinsKa); + int nBinsKa = static_cast(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; @@ -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))) { @@ -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; @@ -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 pvec0 = {_pxPr, _pyPr, _pzPr}; + std::array pvec1 = {_pxKa, _pyKa, _pzKa}; + std::array, 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 pvec0rot = {_pxPr, _pyPr, _pzPr}; + std::array pvec1rot = {_pxKaRot, _pyKaRot, _pzKa}; + std::array, 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 { @@ -542,6 +615,8 @@ struct lambdaAnalysis_pb { } using resoCols = soa::Join; + using resoMCCols = soa::Join; + using resoTracks = aod::ResoTracks; void processData(resoCols::iterator const& collision, resoTracks const& tracks) @@ -549,17 +624,32 @@ struct lambdaAnalysis_pb { // 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(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 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(tracks, tracks, mult); // get MC pT-spectra @@ -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; } @@ -605,10 +695,10 @@ struct lambdaAnalysis_pb { if (!pass1 || !pass2) // If we have both decay products continue; + std::array 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