diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index d6b88601788..4d8d08ba6a9 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -850,7 +850,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(KStarMCMass, kstarMCMass, DECLARE_SOA_DYNAMIC_COLUMN(KStarMCY, kstarMCY, [](float photonMCPx, float photonMCPy, float photonMCPz, float kshortMCPx, float kshortMCPy, float kshortMCPz) -> float { - return RecoDecay::y(std::array{photonMCPx + kshortMCPx, photonMCPy + kshortMCPy, photonMCPz + kshortMCPz}, o2::constants::physics::MassKaonNeutral); + return RecoDecay::y(std::array{photonMCPx + kshortMCPx, photonMCPy + kshortMCPy, photonMCPz + kshortMCPz}, o2::constants::physics::MassK0Star892); }); DECLARE_SOA_DYNAMIC_COLUMN(MCPhi, mcphi, //! Phi in the range [0, 2pi) diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index 73d4997f821..e43f0424ebd 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -1620,13 +1620,13 @@ struct sigma0builder { // Fill with properties GenInfo.IsPrimary = mcParticle.isPhysicalPrimary(); - GenInfo.IsV0Lambda = mcParticle.pdgCode() == PDG_t::kLambda0; // 3122 - GenInfo.IsV0AntiLambda = mcParticle.pdgCode() == PDG_t::kLambda0Bar; //-3122 - GenInfo.IsV0KShort = mcParticle.pdgCode() == PDG_t::kK0Short; // 310 - GenInfo.IsPi0 = mcParticle.pdgCode() == PDG_t::kPi0; // 111; - GenInfo.IsSigma0 = mcParticle.pdgCode() == PDG_t::kSigma0; // PDG_t::kSigma0 - GenInfo.IsAntiSigma0 = mcParticle.pdgCode() == PDG_t::kSigma0Bar; //-3212 - GenInfo.IsKStar = mcParticle.pdgCode() == o2::constants::physics::Pdg::kK0Star892; // 313; + GenInfo.IsV0Lambda = mcParticle.pdgCode() == PDG_t::kLambda0; // 3122 + GenInfo.IsV0AntiLambda = mcParticle.pdgCode() == PDG_t::kLambda0Bar; //-3122 + GenInfo.IsV0KShort = std::abs(mcParticle.pdgCode()) == PDG_t::kK0Short; // 310 + GenInfo.IsPi0 = mcParticle.pdgCode() == PDG_t::kPi0; // 111; + GenInfo.IsSigma0 = mcParticle.pdgCode() == PDG_t::kSigma0; // PDG_t::kSigma0 + GenInfo.IsAntiSigma0 = mcParticle.pdgCode() == PDG_t::kSigma0Bar; //-3212 + GenInfo.IsKStar = std::abs(mcParticle.pdgCode()) == o2::constants::physics::Pdg::kK0Star892; // 313; GenInfo.IsProducedByGenerator = mcParticle.producedByGenerator(); GenInfo.MCProcess = mcParticle.getProcess(); GenInfo.MCPt = mcParticle.pt(); diff --git a/PWGLF/Tasks/Resonances/k892hadronphoton.cxx b/PWGLF/Tasks/Resonances/k892hadronphoton.cxx index fa6def1b51b..62a32bcfc01 100644 --- a/PWGLF/Tasks/Resonances/k892hadronphoton.cxx +++ b/PWGLF/Tasks/Resonances/k892hadronphoton.cxx @@ -809,6 +809,8 @@ struct k892hadronphoton { int trkCode = 10; // 1: TPC-only, 2: TPC+Something, 3: ITS-Only, 4: ITS+TPC + Something, 10: anything else auto photonPosTrackCode = kstar.photonPosTrackCode(); auto photonNegTrackCode = kstar.photonNegTrackCode(); + auto kshortPosTrackCode = kstar.kshortPosTrackCode(); + auto kshortNegTrackCode = kstar.kshortNegTrackCode(); if (isGamma) { if (photonPosTrackCode == DauTPCOnly && photonNegTrackCode == DauTPCOnly) @@ -820,13 +822,13 @@ struct k892hadronphoton { if (photonPosTrackCode == DauITSTrackerOnly || photonNegTrackCode == DauITSTrackerOnly) trkCode = 4; } else { - if (photonPosTrackCode == DauTPCOnly && photonNegTrackCode == DauTPCOnly) + if (kshortPosTrackCode == DauTPCOnly && kshortNegTrackCode == DauTPCOnly) trkCode = 1; - if ((photonPosTrackCode != DauTPCOnly && photonNegTrackCode == DauTPCOnly) || (photonPosTrackCode == DauTPCOnly && photonNegTrackCode != DauTPCOnly)) + if ((kshortPosTrackCode != DauTPCOnly && kshortNegTrackCode == DauTPCOnly) || (kshortPosTrackCode == DauTPCOnly && kshortNegTrackCode != DauTPCOnly)) trkCode = 2; - if (photonPosTrackCode == DauITSTPC && photonNegTrackCode == DauITSTPC) + if (kshortPosTrackCode == DauITSTPC && kshortNegTrackCode == DauITSTPC) trkCode = 3; - if (photonPosTrackCode == DauITSTrackerOnly || photonNegTrackCode == DauITSTrackerOnly) + if (kshortPosTrackCode == DauITSTrackerOnly || kshortNegTrackCode == DauITSTrackerOnly) trkCode = 4; } @@ -899,7 +901,7 @@ struct k892hadronphoton { if (std::abs(kshortPDGCodeGrandMother) == o2::constants::physics::Pdg::kK0Star892 && std::abs(photonPDGCodeGrandMother) == o2::constants::physics::Pdg::kK0Star892 && (photonGlobalIndexGrandMother == kshortGlobalIndexGrandMother)) // K*(892)0 histos.fill(HIST(MainDir[mode]) + HIST("/MC/BkgStudy/h2dPtVsMassKStar_KStarPionKaon"), kstarpT, kstarMass); - if (photonGlobalIndexGrandMother == kshortGlobalIndexGrandMother) { + if ((photonGlobalIndexGrandMother == kshortGlobalIndexGrandMother) && (photonPDGCodeGrandMother != o2::constants::physics::Pdg::kK0Star892) && (photonPDGCodeGrandMother != PDG_t::kProton)) { histos.fill(HIST(MainDir[mode]) + HIST("/MC/BkgStudy/h2dPtVsMassKStar_PionKaon"), kstarpT, kstarMass); histos.fill(HIST(MainDir[mode]) + HIST("/MC/BkgStudy/h2dPionKaonMother"), kshortPDGCodeGrandMother, photonPDGCodeGrandMother); } diff --git a/PWGLF/Tasks/Resonances/k892hadronphotonBkg.cxx b/PWGLF/Tasks/Resonances/k892hadronphotonBkg.cxx index 779757a59c0..239897406f0 100644 --- a/PWGLF/Tasks/Resonances/k892hadronphotonBkg.cxx +++ b/PWGLF/Tasks/Resonances/k892hadronphotonBkg.cxx @@ -62,6 +62,14 @@ static const std::vector photonSels = {"No Sel", "Mass", "Y", "Neg static const std::vector kshortSels = {"No Sel", "Mass", "Y", "Neg Eta", "Pos Eta", "DCAToPV", "Radius", "Z", "DCADau", "Armenteros", "CosPA", "TPCCR", "ITSNCls", "Lifetime", "TPC NSigma"}; +static const std::vector lambdaSels = {"No Sel", "Radius", "Z", "DCADau", "Armenteros", + "CosPA", "Y + Dau Eta", "TPCCR", "ITSNCls", "Lifetime", + "PID", "DCAToPV", "Mass"}; + +enum BkgResonance { + kResoKStar = 0, // K*(892)0 -> K0S + gamma + kResoLambdaStar = 1 // Lambda(1520) -> Lambda + gamma +}; struct k892hadronphotonBkg { Service ccdb{}; @@ -81,6 +89,7 @@ struct k892hadronphotonBkg { Configurable fIRCrashOnNull{"fIRCrashOnNull", false, "Flag to avoid CTP RateFetcher crash."}; Configurable irSource{"irSource", "T0VTX", "Estimator of the interaction rate (Recommended: pp --> T0VTX, Pb-Pb --> ZNC hadronic)"}; + // KStar(892) -> K0S + gamma background struct : ConfigurableGroup { std::string prefix = "kstarBkgConfig"; Configurable doSameEvtRotation{"doSameEvtRotation", false, "Same-event rotational background"}; @@ -95,6 +104,19 @@ struct k892hadronphotonBkg { Configurable rotGamma{"rotGamma", false, "Flag to rotate the photon direction"}; } kstarBkgConfig; + // Lambda(1520) -> Lambda + gamma background + struct : ConfigurableGroup { + std::string prefix = "lambdaStarBkgConfig"; + Configurable doSameEvtRotation{"doSameEvtRotation", false, "Same-event rotational background"}; + Configurable doEvtMixing{"doEvtMixing", false, "Mixed-event background"}; + Configurable lstarMaxOPAngle{"lstarMaxOPAngle", 7.f, "Max opening angle (rad)"}; + Configurable lstarMaxRap{"lstarMaxRap", 0.5f, "Max |y(#Lambda(1520))|"}; + Configurable nBkgRot{"nBkgRot", 3, "Rotations per pair (rotational bkg)"}; + Configurable rotationalCut{"rotationalCut", 10, "theta band: [pi - pi/cut, pi + pi/cut]"}; + Configurable rotationalFactor{"rotationalFactor", 1.f, "Factor to scale the angle of rotation (rotationalFactor * PI)"}; + Configurable rotGamma{"rotGamma", false, "Flag to rotate the photon direction"}; + } lstarBkgConfig; + ConfigurableAxis axisVertexMixBkg{"axisVertexMixBkg", {VARIABLE_WIDTH, -10.f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "z-vertex bins for mixing"}; ConfigurableAxis axisCentralityMixBkg{"axisCentralityMixBkg", {VARIABLE_WIDTH, 0.0f, 1.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 110.0f}, "centrality bins for mixing"}; @@ -181,12 +203,46 @@ struct k892hadronphotonBkg { Configurable kshortMaxTPCNSigmas{"kshortMaxTPCNSigmas", 1e+9, "Max |TPC NSigma| (pion hypothesis) for K0S daughters"}; } kshortSelections; + //// Lambda criteria:: + struct : ConfigurableGroup { + std::string prefix = "lambdaSelections"; // JSON group name + Configurable Lambda_MLThreshold{"Lambda_MLThreshold", 0.1, "Decision Threshold value to select lambdas"}; + Configurable AntiLambda_MLThreshold{"AntiLambda_MLThreshold", 0.1, "Decision Threshold value to select antilambdas"}; + Configurable LambdaMinDCANegToPv{"LambdaMinDCANegToPv", .05, "min DCA Neg To PV (cm)"}; + Configurable LambdaMinDCAPosToPv{"LambdaMinDCAPosToPv", .05, "min DCA Pos To PV (cm)"}; + Configurable ALambdaMinDCANegToPv{"ALambdaMinDCANegToPv", .05, "min DCA Neg To PV (cm)"}; + Configurable ALambdaMinDCAPosToPv{"ALambdaMinDCAPosToPv", .05, "min DCA Pos To PV (cm)"}; + Configurable LambdaMaxDCAV0Dau{"LambdaMaxDCAV0Dau", 2.5, "Max DCA V0 Daughters (cm)"}; + Configurable LambdaMinv0radius{"LambdaMinv0radius", 0.0, "Min V0 radius (cm)"}; + Configurable LambdaMaxv0radius{"LambdaMaxv0radius", 40, "Max V0 radius (cm)"}; + Configurable LambdaMinQt{"LambdaMinQt", 0.01, "Min lambda qt value (AP plot) (GeV/c)"}; + Configurable LambdaMaxQt{"LambdaMaxQt", 0.17, "Max lambda qt value (AP plot) (GeV/c)"}; + Configurable LambdaMinAlpha{"LambdaMinAlpha", 0.25, "Min lambda alpha absolute value (AP plot)"}; + Configurable LambdaMaxAlpha{"LambdaMaxAlpha", 1.0, "Max lambda alpha absolute value (AP plot)"}; + Configurable LambdaMinv0cospa{"LambdaMinv0cospa", 0.95, "Min V0 CosPA"}; + Configurable LambdaMaxLifeTime{"LambdaMaxLifeTime", 30, "Max lifetime"}; + Configurable LambdaWindow{"LambdaWindow", 0.015, "Mass window around expected (in GeV/c2)"}; + Configurable LambdaMinRapidity{"LambdaMinRapidity", -0.5, "v0 min rapidity"}; + Configurable LambdaMaxRapidity{"LambdaMaxRapidity", 0.5, "v0 max rapidity"}; + Configurable LambdaMinDauEta{"LambdaMinDauEta", -0.8, "Min pseudorapidity of daughter tracks"}; + Configurable LambdaMaxDauEta{"LambdaMaxDauEta", 0.8, "Max pseudorapidity of daughter tracks"}; + Configurable LambdaMinZ{"LambdaMinZ", -240, "Min lambda decay point z value (cm)"}; + Configurable LambdaMaxZ{"LambdaMaxZ", 240, "Max lambda decay point z value (cm)"}; + Configurable fselLambdaTPCPID{"fselLambdaTPCPID", true, "Flag to select lambda-like candidates using TPC NSigma."}; + Configurable LambdaMaxTPCNSigmas{"LambdaMaxTPCNSigmas", 1e+9, "Max TPC NSigmas for daughters"}; + Configurable LambdaMinTPCCrossedRows{"LambdaMinTPCCrossedRows", 50, "Min daughter TPC Crossed Rows"}; + Configurable LambdaMinITSclusters{"LambdaMinITSclusters", 1, "minimum ITS clusters"}; + Configurable LambdaRejectPosITSafterburner{"LambdaRejectPosITSafterburner", false, "reject positive track formed out of afterburner ITS tracks"}; + Configurable LambdaRejectNegITSafterburner{"LambdaRejectNegITSafterburner", false, "reject negative track formed out of afterburner ITS tracks"}; + } lambdaSelections; + struct : ConfigurableGroup { // base properties std::string prefix = "axisConfig"; // JSON group name ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "pt axis for analysis"}; ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 110.0f}, "Centrality"}; ConfigurableAxis axisKStarMass{"axisKStarMass", {500, 0.6f, 1.6f}, "M_{K^{*}} (GeV/c^{2})"}; + ConfigurableAxis axisLambdaStarMass{"axisLambdaStarMass", {500, 1.1f, 2.1f}, "M_{#Lambda(1520)} (GeV/c^{2})"}; ConfigurableAxis axisIRBinning{"axisIRBinning", {151, -10, 1500}, "Binning for the interaction rate (kHz)"}; ConfigurableAxis axisCandSel{"axisCandSel", {15, 0.5f, +15.5f}, "Candidate Selection"}; } axisConfig; @@ -244,6 +300,10 @@ struct k892hadronphotonBkg { for (size_t i = 0; i < kshortSels.size(); ++i) histos.get(HIST("KShortSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(i + 1, kshortSels[i].c_str()); + histos.add("LambdaSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisConfig.axisCandSel}); + for (size_t i = 0; i < lambdaSels.size(); ++i) + histos.get(HIST("LambdaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(i + 1, lambdaSels[i].c_str()); + if (kstarBkgConfig.doSameEvtRotation || kstarBkgConfig.doEvtMixing) { histos.add("KStarBkg/hDeltaCollision", "hDeltaCollision", kTH1D, {{2000, -1000.f, 1000.f}}); histos.add("KStarBkg/h2dCentralityCollPair", "h2dCentralityCollPair", kTH2D, {axisConfig.axisCentrality, axisConfig.axisCentrality}); @@ -258,6 +318,22 @@ struct k892hadronphotonBkg { histos.add("KStarBkg/h3dMixedKStarMassVsPt", "h3dMixedKStarMassVsPt", kTH3D, {axisConfig.axisCentrality, axisConfig.axisPt, axisConfig.axisKStarMass}); histos.add("KStarBkg/h3dMixedKStarPtVsOPAngle", "h3dMixedKStarPtVsOPAngle", kTH3D, {{140, 0.f, 7.f}, axisConfig.axisPt, axisConfig.axisKStarMass}); } + + // Lambda(1520) -> Lambda + gamma + if (lstarBkgConfig.doSameEvtRotation || lstarBkgConfig.doEvtMixing) { + histos.add("LambdaStarBkg/hDeltaCollision", "hDeltaCollision", kTH1D, {{2000, -1000.f, 1000.f}}); + histos.add("LambdaStarBkg/h2dCentralityCollPair", "h2dCentralityCollPair", kTH2D, {axisConfig.axisCentrality, axisConfig.axisCentrality}); + } + if (lstarBkgConfig.doSameEvtRotation) { + histos.add("LambdaStarBkg/h2dRotLambdaStarMassVsPt", "h2dRotLambdaStarMassVsPt", kTH2D, {axisConfig.axisLambdaStarMass, axisConfig.axisPt}); + histos.add("LambdaStarBkg/h3dRotLambdaStarMassVsPt", "h3dRotLambdaStarMassVsPt", kTH3D, {axisConfig.axisCentrality, axisConfig.axisPt, axisConfig.axisLambdaStarMass}); + histos.add("LambdaStarBkg/h3dRotLambdaStarPtVsOPAngle", "h3dRotLambdaStarPtVsOPAngle", kTH3D, {{140, 0.f, 7.f}, axisConfig.axisPt, axisConfig.axisLambdaStarMass}); + } + if (lstarBkgConfig.doEvtMixing) { + histos.add("LambdaStarBkg/h2dMixedLambdaStarMassVsPt", "h2dMixedLambdaStarMassVsPt", kTH2D, {axisConfig.axisLambdaStarMass, axisConfig.axisPt}); + histos.add("LambdaStarBkg/h3dMixedLambdaStarMassVsPt", "h3dMixedLambdaStarMassVsPt", kTH3D, {axisConfig.axisCentrality, axisConfig.axisPt, axisConfig.axisLambdaStarMass}); + histos.add("LambdaStarBkg/h3dMixedLambdaStarPtVsOPAngle", "h3dMixedLambdaStarPtVsOPAngle", kTH3D, {{140, 0.f, 7.f}, axisConfig.axisPt, axisConfig.axisLambdaStarMass}); + } } //_______________________________________________ @@ -571,19 +647,137 @@ struct k892hadronphotonBkg { } //_______________________________________________ - // Compute same-event rotational background for K* within a single collision - template + // Process Lambda candidate + template + bool processLambdaCandidate(TV0Object const& lambda, TCollision const& collision) + { + // V0 type selection (matching builder-level) + if (lambda.v0Type() != 1) + return false; + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 1.); + if ((lambda.v0radius() < lambdaSelections.LambdaMinv0radius) || (lambda.v0radius() > lambdaSelections.LambdaMaxv0radius)) + return false; + + // Decay-point z window (builder-level) + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 2.); + if ((lambda.z() < lambdaSelections.LambdaMinZ) || (lambda.z() > lambdaSelections.LambdaMaxZ)) + return false; + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 3.); + if (std::abs(lambda.dcaV0daughters()) > lambdaSelections.LambdaMaxDCAV0Dau) + return false; + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 4.); + if ((lambda.qtarm() < lambdaSelections.LambdaMinQt) || (lambda.qtarm() > lambdaSelections.LambdaMaxQt)) + return false; + + if ((std::abs(lambda.alpha()) < lambdaSelections.LambdaMinAlpha) || (std::abs(lambda.alpha()) > lambdaSelections.LambdaMaxAlpha)) + return false; + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 5.); + if (lambda.v0cosPA() < lambdaSelections.LambdaMinv0cospa) + return false; + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 6.); + if ((lambda.yLambda() < lambdaSelections.LambdaMinRapidity) || (lambda.yLambda() > lambdaSelections.LambdaMaxRapidity)) + return false; + if ((lambda.positiveeta() < lambdaSelections.LambdaMinDauEta) || (lambda.positiveeta() > lambdaSelections.LambdaMaxDauEta)) + return false; + if ((lambda.negativeeta() < lambdaSelections.LambdaMinDauEta) || (lambda.negativeeta() > lambdaSelections.LambdaMaxDauEta)) + return false; + + auto posTrackLambda = lambda.template posTrackExtra_as(); + auto negTrackLambda = lambda.template negTrackExtra_as(); + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 7.); + if ((posTrackLambda.tpcCrossedRows() < lambdaSelections.LambdaMinTPCCrossedRows) || (negTrackLambda.tpcCrossedRows() < lambdaSelections.LambdaMinTPCCrossedRows)) + return false; + + // MinITSCls + reject ITS afterburner tracks if requested + bool posIsFromAfterburner = posTrackLambda.itsChi2PerNcl() < 0; + bool negIsFromAfterburner = negTrackLambda.itsChi2PerNcl() < 0; + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 8.); + if (posTrackLambda.itsNCls() < lambdaSelections.LambdaMinITSclusters && (!lambdaSelections.LambdaRejectPosITSafterburner || posIsFromAfterburner)) + return false; + if (negTrackLambda.itsNCls() < lambdaSelections.LambdaMinITSclusters && (!lambdaSelections.LambdaRejectNegITSafterburner || negIsFromAfterburner)) + return false; + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 9.); + float fLambdaLifeTime = lambda.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassLambda0; + if (fLambdaLifeTime > lambdaSelections.LambdaMaxLifeTime) + return false; + + // Separating lambda and antilambda selections: + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 10.); + if (lambda.alpha() > 0) { // Lambda selection + + // TPC Selection + if (lambdaSelections.fselLambdaTPCPID && (std::abs(posTrackLambda.tpcNSigmaPr()) > lambdaSelections.LambdaMaxTPCNSigmas)) + return false; + if (lambdaSelections.fselLambdaTPCPID && (std::abs(negTrackLambda.tpcNSigmaPi()) > lambdaSelections.LambdaMaxTPCNSigmas)) + return false; + + // DCA Selection + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 11.); + if ((std::abs(lambda.dcapostopv()) < lambdaSelections.LambdaMinDCAPosToPv) || (std::abs(lambda.dcanegtopv()) < lambdaSelections.LambdaMinDCANegToPv)) + return false; + + // Mass Selection + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 12.); + if (std::abs(lambda.mLambda() - o2::constants::physics::MassLambda0) > lambdaSelections.LambdaWindow) + return false; + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 13.); + + } else { // AntiLambda selection + + // TPC Selection + if (lambdaSelections.fselLambdaTPCPID && (std::abs(posTrackLambda.tpcNSigmaPi()) > lambdaSelections.LambdaMaxTPCNSigmas)) + return false; + if (lambdaSelections.fselLambdaTPCPID && (std::abs(negTrackLambda.tpcNSigmaPr()) > lambdaSelections.LambdaMaxTPCNSigmas)) + return false; + + // DCA Selection + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 11.); + if ((std::abs(lambda.dcapostopv()) < lambdaSelections.ALambdaMinDCAPosToPv) || (std::abs(lambda.dcanegtopv()) < lambdaSelections.ALambdaMinDCANegToPv)) + return false; + + // Mass Selection + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 12.); + if (std::abs(lambda.mAntiLambda() - o2::constants::physics::MassLambda0) > lambdaSelections.LambdaWindow) + return false; + + histos.fill(HIST("LambdaSel/hSelectionStatistics"), 13.); + } + + return true; + } + + //_______________________________________________ + // Compute same-event rotational background within a single collision. + template void calculateRotBackground(TCollision const& coll, std::vector const& photonIndices, - std::vector const& kshortIndices, + std::vector const& hadronIndices, TV0s const& fullV0s) { - if (photonIndices.empty() || kshortIndices.empty()) + if (photonIndices.empty() || hadronIndices.empty()) return; - const float centrality = doPPAnalysis ? coll.centFT0M() : coll.centFT0C(); - for (const int& kIdx : kshortIndices) { - const auto& kshort = fullV0s.rawIteratorAt(kIdx); + constexpr float HadronMass = (resonance == kResoKStar) ? o2::constants::physics::MassK0Short : o2::constants::physics::MassLambda0; + constexpr float ResonanceMass = (resonance == kResoKStar) ? o2::constants::physics::MassK0Star892 : o2::constants::physics::MassLambda1520; + + const int nBkgRot = (resonance == kResoKStar) ? kstarBkgConfig.nBkgRot.value : lstarBkgConfig.nBkgRot.value; + const int rotationalCut = (resonance == kResoKStar) ? kstarBkgConfig.rotationalCut.value : lstarBkgConfig.rotationalCut.value; + const float rotationalFactor = (resonance == kResoKStar) ? kstarBkgConfig.rotationalFactor.value : lstarBkgConfig.rotationalFactor.value; + const float maxRap = (resonance == kResoKStar) ? kstarBkgConfig.kstarMaxRap.value : lstarBkgConfig.lstarMaxRap.value; + const bool rotGamma = (resonance == kResoKStar) ? kstarBkgConfig.rotGamma.value : lstarBkgConfig.rotGamma.value; + + const float centrality = getCentralityRun3Bkg(coll); + for (const int& hIdx : hadronIndices) { + const auto& hadron = fullV0s.rawIteratorAt(hIdx); for (const int& pIdx : photonIndices) { const auto& photon = fullV0s.rawIteratorAt(pIdx); @@ -594,41 +788,111 @@ struct k892hadronphotonBkg { photon.phi(), o2::constants::physics::MassGamma); - ROOT::Math::PtEtaPhiMVector pKShort(kshort.pt(), - kshort.eta(), - kshort.phi(), - o2::constants::physics::MassK0Short); + ROOT::Math::PtEtaPhiMVector pHadron(hadron.pt(), + hadron.eta(), + hadron.phi(), + HadronMass); - for (int irot = 0; irot < kstarBkgConfig.nBkgRot; ++irot) { - float theta = rotRng.Uniform(kstarBkgConfig.rotationalFactor * o2::constants::math::PI - o2::constants::math::PI / kstarBkgConfig.rotationalCut, - kstarBkgConfig.rotationalFactor * o2::constants::math::PI + o2::constants::math::PI / kstarBkgConfig.rotationalCut); + for (int irot = 0; irot < nBkgRot; ++irot) { + float theta = rotRng.Uniform(rotationalFactor * o2::constants::math::PI - o2::constants::math::PI / rotationalCut, + rotationalFactor * o2::constants::math::PI + o2::constants::math::PI / rotationalCut); - ROOT::Math::PtEtaPhiMVector kRot(kshort.pt(), kshort.eta(), kshort.phi() + theta, o2::constants::physics::MassK0Short); + ROOT::Math::PtEtaPhiMVector hRot(hadron.pt(), hadron.eta(), hadron.phi() + theta, HadronMass); ROOT::Math::PtEtaPhiMVector gRot(photon.pt(), photon.eta(), photon.phi() + theta, o2::constants::physics::MassGamma); - auto kstar = pGamma + kRot; - if (kstarBkgConfig.rotGamma) { - kstar = gRot + pKShort; + auto reso = pGamma + hRot; + if (rotGamma) { + reso = gRot + pHadron; } - float rapidity = RecoDecay::y(std::array{static_cast(kstar.Px()), - static_cast(kstar.Py()), - static_cast(kstar.Pz())}, - o2::constants::physics::MassK0Star892); - if (std::abs(rapidity) > kstarBkgConfig.kstarMaxRap) + float rapidity = RecoDecay::y(std::array{static_cast(reso.Px()), + static_cast(reso.Py()), + static_cast(reso.Pz())}, + ResonanceMass); + if (std::abs(rapidity) > maxRap) continue; - // Opening angle between photon and rotated K0s (QA only, not used as a cut) - double cosOA = pGamma.Vect().Dot(kRot.Vect()) / (pGamma.P() * kRot.P()); - if (kstarBkgConfig.rotGamma) { - cosOA = gRot.Vect().Dot(pKShort.Vect()) / (gRot.P() * pKShort.P()); + // Opening angle between photon and rotated hadron (QA only, not used as a cut) + double cosOA = pGamma.Vect().Dot(hRot.Vect()) / (pGamma.P() * hRot.P()); + if (rotGamma) { + cosOA = gRot.Vect().Dot(pHadron.Vect()) / (gRot.P() * pHadron.P()); } double openAngle = std::acos(cosOA); - histos.fill(HIST("KStarBkg/h2dRotKStarMassVsPt"), kstar.M(), kstar.Pt()); - histos.fill(HIST("KStarBkg/h3dRotKStarMassVsPt"), centrality, kstar.Pt(), kstar.M()); - histos.fill(HIST("KStarBkg/h3dRotKStarPtVsOPAngle"), openAngle, kstar.Pt(), kstar.M()); + if constexpr (resonance == kResoKStar) { + histos.fill(HIST("KStarBkg/h2dRotKStarMassVsPt"), reso.M(), reso.Pt()); + histos.fill(HIST("KStarBkg/h3dRotKStarMassVsPt"), centrality, reso.Pt(), reso.M()); + histos.fill(HIST("KStarBkg/h3dRotKStarPtVsOPAngle"), openAngle, reso.Pt(), reso.M()); + } else { + histos.fill(HIST("LambdaStarBkg/h2dRotLambdaStarMassVsPt"), reso.M(), reso.Pt()); + histos.fill(HIST("LambdaStarBkg/h3dRotLambdaStarMassVsPt"), centrality, reso.Pt(), reso.M()); + histos.fill(HIST("LambdaStarBkg/h3dRotLambdaStarPtVsOPAngle"), openAngle, reso.Pt(), reso.M()); + } + } + } + } + } + + //_______________________________________________ + // Mixed-event pairing: hadrons and photons come from two different collisions. + // Centrality is taken from the reference collision (the first of the pair) + template + void calculateMixedBackground(TCollision const& refColl, + std::vector const& hadronIndices, + std::vector const& photonIndices, + TV0s const& fullV0s) + { + if (hadronIndices.empty() || photonIndices.empty()) + return; + + constexpr float HadronMass = (resonance == kResoKStar) ? o2::constants::physics::MassK0Short : o2::constants::physics::MassLambda0; + constexpr float ResonanceMass = (resonance == kResoKStar) ? o2::constants::physics::MassK0Star892 : o2::constants::physics::MassLambda1520; + + const float maxOPAngle = (resonance == kResoKStar) ? kstarBkgConfig.kstarMaxOPAngle.value : lstarBkgConfig.lstarMaxOPAngle.value; + const float maxRap = (resonance == kResoKStar) ? kstarBkgConfig.kstarMaxRap.value : lstarBkgConfig.lstarMaxRap.value; + + const float centrality = getCentralityRun3Bkg(refColl); + + for (const int& hIdx : hadronIndices) { + const auto& hadron = fullV0s.rawIteratorAt(hIdx); + float hP = std::hypot(hadron.px(), hadron.py(), hadron.pz()); + ROOT::Math::PxPyPzEVector fourMomHadron( + hadron.px(), hadron.py(), hadron.pz(), + std::sqrt(hP * hP + HadronMass * HadronMass)); + + for (const int& pIdx : photonIndices) { + const auto& photon = fullV0s.rawIteratorAt(pIdx); + float pP = std::hypot(photon.px(), photon.py(), photon.pz()); + ROOT::Math::PxPyPzEVector fourMomPhoton( + photon.px(), photon.py(), photon.pz(), pP); + + auto fourMomReso = fourMomPhoton + fourMomHadron; + + double cosOA = fourMomPhoton.Vect().Dot(fourMomHadron.Vect()) / + (fourMomPhoton.P() * fourMomHadron.P()); + double openAngle = std::acos(cosOA); + float mass = fourMomReso.M(); + float pt = fourMomReso.Pt(); + + float rapidity = RecoDecay::y(std::array{static_cast(fourMomReso.Px()), + static_cast(fourMomReso.Py()), + static_cast(fourMomReso.Pz())}, + ResonanceMass); + + if (openAngle > maxOPAngle) + continue; + if (std::abs(rapidity) > maxRap) + continue; + + if constexpr (resonance == kResoKStar) { + histos.fill(HIST("KStarBkg/h2dMixedKStarMassVsPt"), mass, pt); + histos.fill(HIST("KStarBkg/h3dMixedKStarMassVsPt"), centrality, pt, mass); + histos.fill(HIST("KStarBkg/h3dMixedKStarPtVsOPAngle"), openAngle, pt, mass); + } else { + histos.fill(HIST("LambdaStarBkg/h2dMixedLambdaStarMassVsPt"), mass, pt); + histos.fill(HIST("LambdaStarBkg/h3dMixedLambdaStarMassVsPt"), centrality, pt, mass); + histos.fill(HIST("LambdaStarBkg/h3dMixedLambdaStarPtVsOPAngle"), openAngle, pt, mass); } } } @@ -643,14 +907,18 @@ struct k892hadronphotonBkg { } //_______________________________________________ - // Main: same-event rotation + event mixing for K* background (data only) + // Main: same-event rotation + event mixing for the K* and Lambda(1520) backgrounds using BkgBinningType = ColumnBinningPolicy; template - void calculateKStarBkg(TCollisions const& collisions, TV0s const& fullV0s) + void calculateResonanceBkg(TCollisions const& collisions, TV0s const& fullV0s) { - // Per-collision pools of selected photon and K0s V0 indices + // Lambdas are only selected if a Lambda(1520) background was requested + const bool doLambdaStar = lstarBkgConfig.doSameEvtRotation || lstarBkgConfig.doEvtMixing; + + // Per-collision pools of selected photon, K0s and Lambda V0 indices std::vector> photonPool(collisions.size()); std::vector> kshortPool(collisions.size()); + std::vector> lambdaPool(collisions.size()); // V0 grouping by straCollisionId std::vector> v0grouped(collisions.size()); @@ -672,19 +940,28 @@ struct k892hadronphotonBkg { if (processKShortCandidate(v0, coll)) kshortPool[coll.globalIndex()].push_back(v0.globalIndex()); + + if (doLambdaStar && processLambdaCandidate(v0, coll)) + lambdaPool[coll.globalIndex()].push_back(v0.globalIndex()); } // Same-event rotational background if (kstarBkgConfig.doSameEvtRotation) { - calculateRotBackground(coll, - photonPool[coll.globalIndex()], - kshortPool[coll.globalIndex()], - fullV0s); + calculateRotBackground(coll, + photonPool[coll.globalIndex()], + kshortPool[coll.globalIndex()], + fullV0s); + } + if (lstarBkgConfig.doSameEvtRotation) { + calculateRotBackground(coll, + photonPool[coll.globalIndex()], + lambdaPool[coll.globalIndex()], + fullV0s); } } // Event Mixing - if (!kstarBkgConfig.doEvtMixing) + if (!kstarBkgConfig.doEvtMixing && !lstarBkgConfig.doEvtMixing) return; // Build the mixing binning locally: a struct member initialized from a @@ -696,125 +973,49 @@ struct k892hadronphotonBkg { if (coll1.globalIndex() == coll2.globalIndex()) continue; - histos.fill(HIST("KStarBkg/hDeltaCollision"), - coll1.globalIndex() - coll2.globalIndex()); - histos.fill(HIST("KStarBkg/h2dCentralityCollPair"), - getCentralityRun3Bkg(coll1), getCentralityRun3Bkg(coll2)); + if (kstarBkgConfig.doEvtMixing) { + histos.fill(HIST("KStarBkg/hDeltaCollision"), + coll1.globalIndex() - coll2.globalIndex()); + histos.fill(HIST("KStarBkg/h2dCentralityCollPair"), + getCentralityRun3Bkg(coll1), getCentralityRun3Bkg(coll2)); + } + if (lstarBkgConfig.doEvtMixing) { + histos.fill(HIST("LambdaStarBkg/hDeltaCollision"), + coll1.globalIndex() - coll2.globalIndex()); + histos.fill(HIST("LambdaStarBkg/h2dCentralityCollPair"), + getCentralityRun3Bkg(coll1), getCentralityRun3Bkg(coll2)); + } if (std::abs(static_cast(coll1.globalIndex()) - static_cast(coll2.globalIndex())) < kstarBkgConfig.deltaCollision) continue; auto const& photons1 = photonPool[coll1.globalIndex()]; - auto const& kshorts1 = kshortPool[coll1.globalIndex()]; auto const& photons2 = photonPool[coll2.globalIndex()]; - auto const& kshorts2 = kshortPool[coll2.globalIndex()]; - - // K0s(coll1) × γ(coll2) - if (!kshorts1.empty() && !photons2.empty()) { - for (const int& kIdx : kshorts1) { - const auto& kshort = fullV0s.rawIteratorAt(kIdx); - float kP = std::hypot(kshort.px(), kshort.py(), kshort.pz()); - ROOT::Math::PxPyPzEVector fourMomKShort( - kshort.px(), kshort.py(), kshort.pz(), - std::sqrt(kP * kP + - o2::constants::physics::MassK0Short * - o2::constants::physics::MassK0Short)); - - for (const int& pIdx : photons2) { - const auto& photon = fullV0s.rawIteratorAt(pIdx); - float pP = std::hypot(photon.px(), photon.py(), photon.pz()); - ROOT::Math::PxPyPzEVector fourMomPhoton( - photon.px(), photon.py(), photon.pz(), pP); - - auto fourMomKStar = fourMomPhoton + fourMomKShort; - - double cosOA = fourMomPhoton.Vect().Dot(fourMomKShort.Vect()) / - (fourMomPhoton.P() * fourMomKShort.P()); - double openAngle = std::acos(cosOA); - float mass = fourMomKStar.M(); - float pt = fourMomKStar.Pt(); - // Rapidity computed under the K*(892) mass hypothesis (NOT the actual pair - // invariant mass) to match the same-event rotational background and the - // buildKStar signal selection, so the rapidity acceptance is identical for - // signal and all backgrounds. - float rapidity = RecoDecay::y(std::array{static_cast(fourMomKStar.Px()), - static_cast(fourMomKStar.Py()), - static_cast(fourMomKStar.Pz())}, - o2::constants::physics::MassK0Star892); - - if (openAngle > kstarBkgConfig.kstarMaxOPAngle) - continue; - if (std::abs(rapidity) > kstarBkgConfig.kstarMaxRap) - continue; - - histos.fill(HIST("KStarBkg/h2dMixedKStarMassVsPt"), mass, pt); - histos.fill(HIST("KStarBkg/h3dMixedKStarMassVsPt"), - getCentralityRun3Bkg(coll1), pt, mass); - histos.fill(HIST("KStarBkg/h3dMixedKStarPtVsOPAngle"), - openAngle, pt, mass); - } - } + + if (kstarBkgConfig.doEvtMixing) { + // K0s(coll1) × γ(coll2) and γ(coll1) × K0s(coll2) + calculateMixedBackground(coll1, kshortPool[coll1.globalIndex()], photons2, fullV0s); + calculateMixedBackground(coll1, kshortPool[coll2.globalIndex()], photons1, fullV0s); } - // γ(coll1) × K0s(coll2) - if (!photons1.empty() && !kshorts2.empty()) { - for (const int& pIdx : photons1) { - const auto& photon = fullV0s.rawIteratorAt(pIdx); - float pP = std::hypot(photon.px(), photon.py(), photon.pz()); - ROOT::Math::PxPyPzEVector fourMomPhoton( - photon.px(), photon.py(), photon.pz(), pP); - - for (const int& kIdx : kshorts2) { - const auto& kshort = fullV0s.rawIteratorAt(kIdx); - float kP = std::hypot(kshort.px(), kshort.py(), kshort.pz()); - ROOT::Math::PxPyPzEVector fourMomKShort( - kshort.px(), kshort.py(), kshort.pz(), - std::sqrt(kP * kP + - o2::constants::physics::MassK0Short * - o2::constants::physics::MassK0Short)); - - auto fourMomKStar = fourMomPhoton + fourMomKShort; - - double cosOA = fourMomPhoton.Vect().Dot(fourMomKShort.Vect()) / - (fourMomPhoton.P() * fourMomKShort.P()); - double openAngle = std::acos(cosOA); - float mass = fourMomKStar.M(); - float pt = fourMomKStar.Pt(); - // Rapidity computed under the K*(892) mass hypothesis (NOT the actual pair - // invariant mass) to match the same-event rotational background and the - // buildKStar signal selection, so the rapidity acceptance is identical for - // signal and all backgrounds. - float rapidity = RecoDecay::y(std::array{static_cast(fourMomKStar.Px()), - static_cast(fourMomKStar.Py()), - static_cast(fourMomKStar.Pz())}, - o2::constants::physics::MassK0Star892); - - if (openAngle > kstarBkgConfig.kstarMaxOPAngle) - continue; - if (std::abs(rapidity) > kstarBkgConfig.kstarMaxRap) - continue; - - histos.fill(HIST("KStarBkg/h2dMixedKStarMassVsPt"), mass, pt); - histos.fill(HIST("KStarBkg/h3dMixedKStarMassVsPt"), - getCentralityRun3Bkg(coll1), pt, mass); - histos.fill(HIST("KStarBkg/h3dMixedKStarPtVsOPAngle"), - openAngle, pt, mass); - } - } + if (lstarBkgConfig.doEvtMixing) { + // Λ(coll1) × γ(coll2) and γ(coll1) × Λ(coll2) + calculateMixedBackground(coll1, lambdaPool[coll1.globalIndex()], photons2, fullV0s); + calculateMixedBackground(coll1, lambdaPool[coll2.globalIndex()], photons1, fullV0s); } } } //_______________________________________________ - // Data process: same-event rotational + mixed-event K* background + // Data process: same-event rotational + mixed-event K* and Lambda(1520) background void processKStarBkg(soa::Join const& collisions, V0StandardDerivedDatas const& fullV0s, dauTracks const&) { - calculateKStarBkg(collisions, fullV0s); + calculateResonanceBkg(collisions, fullV0s); } - PROCESS_SWITCH(k892hadronphotonBkg, processKStarBkg, "Compute K* same-event rotational and mixed-event background (data)", true); + PROCESS_SWITCH(k892hadronphotonBkg, processKStarBkg, "Compute K* and Lambda(1520) same-event rotational and mixed-event background (data)", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)