From b6c8e03f3d216fca029801523118bae77e7c2e40 Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Sat, 16 May 2026 21:41:23 +0200 Subject: [PATCH 1/9] Add histograms to compute event loss and event splitting --- .../Strangeness/strangenessInJetsIons.cxx | 29 ++++++++++++------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index ffb3b001f4d..dfeb031cccd 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -346,9 +346,11 @@ struct StrangenessInJetsIons { // Add histogram to store multiplicity of the event registryMC.add("number_of_events_vsmultiplicity_gen", "number of events vs multiplicity", HistType::kTH1D, {{101, 0, 101, "Multiplicity percentile"}}); + registryMC.add("number_of_events_vsmultiplicity_gen_w_reco", "number of events vs multiplicity (gen with reco)", HistType::kTH1D, {{101, 0, 101, "Multiplicity percentile"}}); // For MB registryMC.add("number_of_events_vsmultiplicity_gen_MB", "number of events vs multiplicity (MB)", HistType::kTH1D, {{101, 0, 101, "Multiplicity percentile"}}); + registryMC.add("number_of_events_vsmultiplicity_gen_w_reco_MB", "number of events vs multiplicity (MB gen with reco)", HistType::kTH1D, {{101, 0, 101, "Multiplicity percentile"}}); // Jet counters registryMC.add("n_jets_vs_mult_pT_mc_gen", "n_jets_vs_mult_pT_mc_gen", HistType::kTH2F, {multAxis, ptJetAxis}); @@ -1187,9 +1189,9 @@ struct StrangenessInJetsIons { return std::abs(yParticle) < rapidityCut; } - void FillFullEventHistoMCGEN(aod::McParticle const& particle, - const double& genMultiplicity, - bool hasReco = false) + void FillMBEventHistoMCGEN(aod::McParticle const& particle, + const double& genMultiplicity, + bool hasReco = false) { if (particle.isPhysicalPrimary() && passedRapidityCut(particle.y(), configV0.rapidityMax)) { switch (particle.pdgCode()) { @@ -1275,12 +1277,12 @@ struct StrangenessInJetsIons { // Fill Minimum Bias histograms template - void FillFullEventHistoMCREC(MCRecoCollision collision, - aod::McParticles const& mcParticles, - V0PerColl const& v0sPerColl, - CascPerColl const& cascPerColl, - TracksPerColl const& tracksPerColl, - const double& multiplicity) + void FillMBEventHistoMCREC(MCRecoCollision collision, + aod::McParticles const& mcParticles, + V0PerColl const& v0sPerColl, + CascPerColl const& cascPerColl, + TracksPerColl const& tracksPerColl, + const double& multiplicity) { // V0 particles if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) { @@ -1913,6 +1915,8 @@ struct StrangenessInJetsIons { genMultiplicity = collision.centFT0M(); } registryMC.fill(HIST("number_of_events_vsmultiplicity_gen_MB"), genMultiplicity); + if (hasReco) + registryMC.fill(HIST("number_of_events_vsmultiplicity_gen_w_reco_MB"), genMultiplicity); // MC particles per collision auto mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, collision.globalIndex()); @@ -1930,7 +1934,7 @@ struct StrangenessInJetsIons { if (particle.eta() < configTracks.etaMin || particle.eta() > configTracks.etaMax || particle.pt() < minPtParticle) continue; - FillFullEventHistoMCGEN(particle, genMultiplicity, hasReco); + FillMBEventHistoMCGEN(particle, genMultiplicity, hasReco); // Select physical primary particles or HF decay products if (!isPhysicalPrimaryOrFromHF(particle, mcParticles)) @@ -1948,6 +1952,7 @@ struct StrangenessInJetsIons { if (fjParticles.size() < 1) continue; registryMC.fill(HIST("number_of_events_mc_gen"), 2.5); + // Cluster MC particles into jets using anti-kt algorithm fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); @@ -1974,6 +1979,8 @@ struct StrangenessInJetsIons { countSelJet++; registryMC.fill(HIST("number_of_events_mc_gen"), 3.5); registryMC.fill(HIST("number_of_events_vsmultiplicity_gen"), genMultiplicity); + if (hasReco) + registryMC.fill(HIST("number_of_events_vsmultiplicity_gen_w_reco"), genMultiplicity); // Fill jet counter registryMC.fill(HIST("n_jets_vs_mult_pT_mc_gen"), genMultiplicity, jetMinusBkg.pt()); @@ -2229,7 +2236,7 @@ struct StrangenessInJetsIons { auto tracksPerColl = mcTracks.sliceBy(perCollisionTrk, collision.globalIndex()); const auto& mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, mcCollision.globalIndex()); - FillFullEventHistoMCREC(collision, mcParticles, v0sPerColl, + FillMBEventHistoMCREC(collision, mcParticles, v0sPerColl, cascPerColl, tracksPerColl, multiplicity); // Loop over reconstructed tracks From 69ad2817d650757420b65e5a37bb4f5edff1d6f2 Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Sat, 16 May 2026 22:26:07 +0200 Subject: [PATCH 2/9] Add process function for MC model predictions --- .../Strangeness/strangenessInJetsIons.cxx | 293 +++++++++++++++++- 1 file changed, 290 insertions(+), 3 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index dfeb031cccd..f4a4bc55173 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -334,7 +334,10 @@ struct StrangenessInJetsIons { } // Histograms for mc generated - if (doprocessMCgenerated) { + if (doprocessMCgenerated || doprocessMCmodels) { + if (doprocessMCgenerated && doprocessMCmodels) { + LOG(fatal) << "processMCgenerated and processMCmodels cannot be selected at the same time. Exit." << endl; + } // Event counter registryMC.add("number_of_events_mc_gen", "number of gen events in mc", HistType::kTH1D, {{10, 0, 10, "Event Cuts"}}); @@ -1952,7 +1955,6 @@ struct StrangenessInJetsIons { if (fjParticles.size() < 1) continue; registryMC.fill(HIST("number_of_events_mc_gen"), 2.5); - // Cluster MC particles into jets using anti-kt algorithm fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); @@ -2237,7 +2239,7 @@ struct StrangenessInJetsIons { const auto& mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, mcCollision.globalIndex()); FillMBEventHistoMCREC(collision, mcParticles, v0sPerColl, - cascPerColl, tracksPerColl, multiplicity); + cascPerColl, tracksPerColl, multiplicity); // Loop over reconstructed tracks for (auto const& track : tracksPerColl) { @@ -2570,6 +2572,291 @@ struct StrangenessInJetsIons { } } PROCESS_SWITCH(StrangenessInJetsIons, processMCreconstructed, "process reconstructed events", false); + + void processMCmodels(soa::Join const& collisions, + aod::McParticles const& mcParticles) + { + // Define per-event particle containers + std::vector fjParticles; + std::vector strHadronMomentum; + std::vector pdg; + + // Jet and area definitions + fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet); + fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0)); + + // Loop over MC collisions + for (const auto& collision : collisions) { + // Clear containers at the start of the event loop + fjParticles.clear(); + strHadronMomentum.clear(); + pdg.clear(); + + // Fill event counter before any selection + registryMC.fill(HIST("number_of_events_mc_gen"), 0.5); + + // Require vertex position within the allowed z range + // if (std::fabs(collision.posZ()) > zVtx) + // continue; + + // Fill event counter after selection on z-vertex + registryMC.fill(HIST("number_of_events_mc_gen"), 1.5); + + // Multiplicity of generated event retrived from corresponding MC RECO collision + // float genMultiplicity = multiplicity; + + // Multiplicity of generated event + float genMultiplicity; + if (centrEstimator == 0) { + genMultiplicity = collision.centFT0C(); + } else { + genMultiplicity = collision.centFT0M(); + } + registryMC.fill(HIST("number_of_events_vsmultiplicity_gen_MB"), genMultiplicity); + + // MC particles per collision + auto mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, collision.globalIndex()); + + // Loop over all MC particles and select physical primaries within acceptance + for (const auto& particle : mcParticlesPerColl) { + // Store properties of strange hadrons + int pdgAbs = std::abs(particle.pdgCode()); + if (particle.isPhysicalPrimary() && (pdgAbs == kK0Short || pdgAbs == kLambda0 || pdgAbs == kXiMinus || pdgAbs == kOmegaMinus)) { // TODO: add protons, kaons, pions + pdg.emplace_back(particle.pdgCode()); + strHadronMomentum.emplace_back(particle.px(), particle.py(), particle.pz()); + } + + double minPtParticle = 0.1f; + if (particle.eta() < configTracks.etaMin || particle.eta() > configTracks.etaMax || particle.pt() < minPtParticle) + continue; + + // false --> do not fill histograms with reco events + FillMBEventHistoMCGEN(particle, genMultiplicity, false); + + // Select physical primary particles or HF decay products + if (!isPhysicalPrimaryOrFromHF(particle, mcParticles)) + continue; + + // Build 4-momentum assuming charged pion mass + static constexpr float kMassPionChargedSquared = o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged; + const double energy = std::sqrt(particle.p() * particle.p() + kMassPionChargedSquared); + fastjet::PseudoJet fourMomentum(particle.px(), particle.py(), particle.pz(), energy); + fourMomentum.set_user_index(particle.pdgCode()); + fjParticles.emplace_back(fourMomentum); + } + + // Skip events with no particles + if (fjParticles.size() < 1) + continue; + registryMC.fill(HIST("number_of_events_mc_gen"), 2.5); + + // Cluster MC particles into jets using anti-kt algorithm + fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); + std::vector jets = fastjet::sorted_by_pt(cs.inclusive_jets()); + + // Estimate background energy density (rho) in perpendicular cone + auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet); + + // Loop over clustered jets + int countSelJet = 0; // number of selected jets + for (const auto& jet : jets) { + + // Jet must be fully contained in acceptance + if ((std::fabs(jet.eta()) + rJet) > (configTracks.etaMax - deltaEtaEdge)) + continue; + + // Subtract background energy from jet + auto jetForSub = jet; + fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp); + + // Apply jet pT threshold + if (jetMinusBkg.pt() < minJetPt) + continue; + countSelJet++; + registryMC.fill(HIST("number_of_events_mc_gen"), 3.5); + registryMC.fill(HIST("number_of_events_vsmultiplicity_gen"), genMultiplicity); + + // Fill jet counter + registryMC.fill(HIST("n_jets_vs_mult_pT_mc_gen"), genMultiplicity, jetMinusBkg.pt()); + + // Set up two perpendicular cone axes for underlying event estimation + TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); + double coneRadius = std::sqrt(jet.area() / PI); + TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0); + getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2); + if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) { + continue; + } + + // Loop over strange hadrons + int index = -1; + for (const auto& hadron : strHadronMomentum) { + // Particle index + index++; + + // Compute distance of particles from jet and UE axes + double deltaEtaJet = hadron.Eta() - jetAxis.Eta(); + double deltaPhiJet = getDeltaPhi(hadron.Phi(), jetAxis.Phi()); + double deltaRJet = std::sqrt(deltaEtaJet * deltaEtaJet + deltaPhiJet * deltaPhiJet); + double deltaEtaUe1 = hadron.Eta() - ueAxis1.Eta(); + double deltaPhiUe1 = getDeltaPhi(hadron.Phi(), ueAxis1.Phi()); + double deltaRUe1 = std::sqrt(deltaEtaUe1 * deltaEtaUe1 + deltaPhiUe1 * deltaPhiUe1); + double deltaEtaUe2 = hadron.Eta() - ueAxis2.Eta(); + double deltaPhiUe2 = getDeltaPhi(hadron.Phi(), ueAxis2.Phi()); + double deltaRUe2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2); + + // Select particles inside jet + if (deltaRJet < coneRadius) { + switch (pdg[index]) { + case kK0Short: + if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + registryMC.fill(HIST("K0s_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kLambda0: + if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + registryMC.fill(HIST("Lambda_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kLambda0Bar: + if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + registryMC.fill(HIST("AntiLambda_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kXiMinus: + if (particleOfInterestDict[ParticleOfInterest::kCascades]) { + registryMC.fill(HIST("XiNeg_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kXiPlusBar: + if (particleOfInterestDict[ParticleOfInterest::kCascades]) { + registryMC.fill(HIST("XiPos_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kOmegaMinus: + if (particleOfInterestDict[ParticleOfInterest::kCascades]) { + registryMC.fill(HIST("OmegaNeg_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kOmegaPlusBar: + if (particleOfInterestDict[ParticleOfInterest::kCascades]) { + registryMC.fill(HIST("OmegaPos_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kPiPlus: + if (particleOfInterestDict[ParticleOfInterest::kPions]) { + registryMC.fill(HIST("PionPos_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kKPlus: + if (particleOfInterestDict[ParticleOfInterest::kKaons]) { + registryMC.fill(HIST("KaonPos_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kProton: + if (particleOfInterestDict[ParticleOfInterest::kProtons]) { + registryMC.fill(HIST("ProtonPos_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kPiMinus: + if (particleOfInterestDict[ParticleOfInterest::kPions]) { + registryMC.fill(HIST("PionNeg_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kKMinus: + if (particleOfInterestDict[ParticleOfInterest::kKaons]) { + registryMC.fill(HIST("KaonNeg_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + case kProtonBar: + if (particleOfInterestDict[ParticleOfInterest::kProtons]) { + registryMC.fill(HIST("ProtonNeg_generated_jet"), genMultiplicity, hadron.Pt()); + } + break; + default: + break; + } + } + + // Select particles inside UE cones + if (deltaRUe1 < coneRadius || deltaRUe2 < coneRadius) { + switch (pdg[index]) { + case kK0Short: + if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + registryMC.fill(HIST("K0s_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kLambda0: + if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + registryMC.fill(HIST("Lambda_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kLambda0Bar: + if (particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + registryMC.fill(HIST("AntiLambda_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kXiMinus: + if (particleOfInterestDict[ParticleOfInterest::kCascades]) { + registryMC.fill(HIST("XiNeg_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kXiPlusBar: + if (particleOfInterestDict[ParticleOfInterest::kCascades]) { + registryMC.fill(HIST("XiPos_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kOmegaMinus: + if (particleOfInterestDict[ParticleOfInterest::kCascades]) { + registryMC.fill(HIST("OmegaNeg_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kOmegaPlusBar: + if (particleOfInterestDict[ParticleOfInterest::kCascades]) { + registryMC.fill(HIST("OmegaPos_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kPiPlus: + if (particleOfInterestDict[ParticleOfInterest::kPions]) { + registryMC.fill(HIST("PionPos_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kKPlus: + if (particleOfInterestDict[ParticleOfInterest::kKaons]) { + registryMC.fill(HIST("KaonPos_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kProton: + if (particleOfInterestDict[ParticleOfInterest::kProtons]) { + registryMC.fill(HIST("ProtonPos_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kPiMinus: + if (particleOfInterestDict[ParticleOfInterest::kPions]) { + registryMC.fill(HIST("PionNeg_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kKMinus: + if (particleOfInterestDict[ParticleOfInterest::kKaons]) { + registryMC.fill(HIST("KaonNeg_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + case kProtonBar: + if (particleOfInterestDict[ParticleOfInterest::kProtons]) { + registryMC.fill(HIST("ProtonNeg_generated_ue"), genMultiplicity, hadron.Pt()); + } + break; + default: + break; + } + } + } + } + // Fill jet counter + registryMC.fill(HIST("n_jets_vs_mult_mc_gen"), genMultiplicity, countSelJet); + } + } + PROCESS_SWITCH(StrangenessInJetsIons, processMCmodels, "Process MC generated events (with models)", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 7646f277d0b278b3cef4bfb694b9bf5d94bba38b Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Sun, 17 May 2026 23:31:32 +0200 Subject: [PATCH 3/9] Add V0s in jet reconstruction in MC RECO with configurable useV0inJetRec --- .../Strangeness/strangenessInJetsIons.cxx | 173 +++++++++++++++++- 1 file changed, 170 insertions(+), 3 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index f4a4bc55173..10db23367aa 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -22,6 +22,7 @@ #include "PWGJE/Core/JetBkgSubUtils.h" #include "PWGJE/Core/JetUtilities.h" +// #include "PWGJE/Core/JetV0Utilities.h" #include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGLF/DataModel/mcCentrality.h" @@ -68,6 +69,7 @@ #include #include #include +#include #include using namespace std; @@ -130,6 +132,7 @@ struct StrangenessInJetsIons { Configurable triggerName{"triggerName", "fOmega", "Software trigger name"}; Configurable centrEstimator{"centrEstimator", 1, "Select centrality estimator. Options: 0 = FT0C, 1 = FT0M. CCDB objects available only for FT0M."}; Configurable calculateFeeddownMatrix{"calculateFeeddownMatrix", true, "Fill feeddown matrix for Lambda if MC"}; + Configurable useV0inJetRec{"useV0inJetRec", true, "Include V0s in jet reconstruction"}; // Event selection Configurable requireNoSameBunchPileup{"requireNoSameBunchPileup", true, "Require kNoSameBunchPileup selection"}; @@ -563,6 +566,12 @@ struct StrangenessInJetsIons { return deltaPhi; } + // Compute energy + double GetEnergy(double px, double py, double pz, double mass) + { + return std::sqrt(px * px + py * py + pz * pz + mass * mass); + } + // Check if particle is a physical primary or a decay product of a heavy-flavor hadron bool isPhysicalPrimaryOrFromHF(aod::McParticle const& particle, aod::McParticles const& mcParticles) { @@ -1529,6 +1538,145 @@ struct StrangenessInJetsIons { } } + /** + * @brief Identifies the index of the common mother for two MC particles. + * + * @tparam T Type of the MC particle objects + * @param posParticle MC particle associated with the positive daughter + * @param negParticle MC particle associated with the negative daughter + * @return int The index of the common mother, or -1 if no common mother exists. + */ + template + int GetCommonMotherId(aod::McParticles const& mcParticles, + const T& posParticle, const T& negParticle) + { + int motherIdPos = posParticle.mothersIds()[0]; + int motherIdNeg = negParticle.mothersIds()[0]; + + const auto& motherPos = mcParticles.iteratorAt(motherIdPos); + const auto& motherNeg = mcParticles.iteratorAt(motherIdNeg); + if (motherPos != motherNeg) { + return -1; // Return -1 if they originate from different parents + } + + return motherIdPos; // Return mother ID + } + + /** + * Returns true if the track is a daughter of the V0 candidate + * Adapted from: PWGJE/Core/JetV0Utilities.h + * + * @param track track that is being checked + * @param candidate V0 candidate that is being checked + */ + template + bool isV0DaughterTrack(T& track, U& candidate) + { + if (candidate.posTrackId() == track.globalIndex() || candidate.negTrackId() == track.globalIndex()) { + return true; + } else { + return false; + } + } + + /** + * @brief Add V0s as input for the jet finder algorithm in MC reco (detector level) + * + * @tparam Track + * @tparam V0PerColl The container type holding the sliced V0 candidates for the current collision. + * + * @param[in,out] fjInput Vector of FastJet PseudoJets where valid V0s will be appended. + * @param[in] fjTracks Vector containing tracks already selected for jet finder input. + * @param[in] v0sPerColl Sliced V0 candidates belonging to this specific collision. + * @param[in] mcParticles Global MC particles table used to look up mother particle truth information. + * @param[in] vtxPos TVector3 object representing the vertex position. + */ + template + void AddV0sForJetReconstructionMC(std::vector& fjInput, + const std::vector& fjTracks, + V0PerColl const& v0sPerColl, + aod::McParticles const& mcParticles, + const TVector3& vtxPos) + { + // Vector of labels: if true remove the element from fjInput + std::vector isTrackReplaced(fjTracks.size(), false); + std::vector v0PseudoJets; // V0s to use as input for jet finder + + for (const auto& v0 : v0sPerColl) { + const auto& pos = v0.template posTrack_as(); + const auto& neg = v0.template negTrack_as(); + + // Get MC particles + if (!pos.has_mcParticle() || !neg.has_mcParticle()) { + continue; + } + const auto& posParticle = pos.template mcParticle_as(); + const auto& negParticle = neg.template mcParticle_as(); + if (!posParticle.has_mothers() || !negParticle.has_mothers()) { + continue; + } + + // Select particles originating from the same parent + int motherId = GetCommonMotherId(mcParticles, posParticle, negParticle); + if (motherId < 0) + continue; // if motherId is -1, it means no common mother was found + const auto& mother = mcParticles.iteratorAt(motherId); + // const bool isPhysPrim = mother.isPhysicalPrimary(); + int pdgParent = mother.pdgCode(); + + /// NOTE: V0s are not required to be primary particles; + bool isK0S = false, isLambda = false, isAntiLambda = false; + // K0s + if (passedK0ShortSelection(v0, pos, neg, vtxPos) && pdgParent == kK0Short) { + // LOG(info) << "[AddV0sForJetReconstructionMC] Add K0S as input for jet finder."; + double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassK0Short); + fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); + v0PseudoJets.emplace_back(fourMomentum); + isK0S = true; + } + // Lambda + if (passedLambdaSelection(v0, pos, neg, vtxPos) && pdgParent == kLambda0) { + // LOG(info) << "[AddV0sForJetReconstructionMC] Add Lambda as input for jet finder."; + double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassLambda0); + fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); + v0PseudoJets.emplace_back(fourMomentum); + isLambda = true; + } + // AntiLambda + if (passedAntiLambdaSelection(v0, pos, neg, vtxPos) && pdgParent == kLambda0Bar) { + // LOG(info) << "[AddV0sForJetReconstructionMC] Add AntiLambda as input for jet finder."; + double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassLambda0Bar); + fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); + v0PseudoJets.emplace_back(fourMomentum); + isAntiLambda = true; + } + + // Exclude daughter tracks in case a V0 is found + bool isV0 = isK0S || isLambda || isAntiLambda; + if (!isV0) + continue; + for (int i = 0; i < fjTracks.size(); ++i) { + if (isV0DaughterTrack(fjTracks[i], v0)) { + // LOG(info) << "[AddV0sForJetReconstructionMC] V0 daughter track found in fjTracks."; + isTrackReplaced[i] = true; + } + } + } + + std::vector cleanFjInput; + cleanFjInput.reserve(fjInput.size()); + for (size_t i = 0; i < fjInput.size(); ++i) { + if (!isTrackReplaced[i]) + cleanFjInput.push_back(fjInput[i]); + } + for (const auto& v0pj : v0PseudoJets) { + cleanFjInput.push_back(v0pj); + } + // LOG(info) << "[AddV0sForJetReconstructionMC] Size fjInput: " << fjInput.size(); + fjInput = std::move(cleanFjInput); + // LOG(info) << "[AddV0sForJetReconstructionMC] Size fjInput dopo move: " << fjInput.size(); + } + // Process data void processData(SelCollisions::iterator const& collision, aod::V0Datas const& fullV0s, aod::CascDataExt const& Cascades, DaughterTracks const& tracks, @@ -2190,6 +2338,9 @@ struct StrangenessInJetsIons { const auto& mcCollision = collision.mcCollision_as>(); + // Vertex position vector + TVector3 vtxPos(collision.posX(), collision.posY(), collision.posZ()); + // Clear containers at the start of the event loop fjParticles.clear(); selectedJet.clear(); @@ -2242,6 +2393,7 @@ struct StrangenessInJetsIons { cascPerColl, tracksPerColl, multiplicity); // Loop over reconstructed tracks + std::vector> fjTracks; for (auto const& track : tracksPerColl) { if (!passedTrackSelectionForJetReconstruction(track)) continue; @@ -2249,7 +2401,14 @@ struct StrangenessInJetsIons { // 4-momentum representation of a particle fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(o2::constants::physics::MassPionCharged)); fjParticles.emplace_back(fourMomentum); + fjTracks.push_back(track); + } + + // Include V0s as tracks for jet reconstruction + if (useV0inJetRec && particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + AddV0sForJetReconstructionMC(fjParticles, fjTracks, v0sPerColl, mcParticles, vtxPos); } + fjTracks.clear(); // Reject empty events if (fjParticles.size() < 1) @@ -2405,6 +2564,17 @@ struct StrangenessInJetsIons { if (pdgParent == 0) continue; + // --- alternative method to avoid double for loop + int motherId = GetCommonMotherId(mcParticles, posParticle, negParticle); + if (motherId < 0) + continue; // if motherId is -1, it means no common mother was found + const auto& motherTest = mcParticles.iteratorAt(motherId); + // const bool isPhysPrimTest = motherTest.isPhysicalPrimary(); + int pdgParentTest = motherTest.pdgCode(); + if (pdgParent != pdgParentTest) + LOG(warning) << "pdgParent != pdgParentTest" << endl; + // --- + // Compute distance from jet and UE axes double deltaEtaJet = v0dir.Eta() - selectedJet[i].Eta(); double deltaPhiJet = getDeltaPhi(v0dir.Phi(), selectedJet[i].Phi()); @@ -2416,9 +2586,6 @@ struct StrangenessInJetsIons { double deltaPhiUe2 = getDeltaPhi(v0dir.Phi(), ue2[i].Phi()); double deltaRue2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2); - // Vertex position vector - TVector3 vtxPos(collision.posX(), collision.posY(), collision.posZ()); - // Fill Armenteros-Podolanski TH2 // registryQC.fill(HIST("ArmenterosPreSel_REC"), v0.alpha(), v0.qtarm()); From c05afac98944d0cd8e17117c7fab0079efe3820d Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Mon, 18 May 2026 09:41:56 +0200 Subject: [PATCH 4/9] Add V0s in jet reconstruction in DATA with configurable useV0inJetRec --- .../Strangeness/strangenessInJetsIons.cxx | 115 +++++++++++++++--- 1 file changed, 100 insertions(+), 15 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index 10db23367aa..72869b1f5e2 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -1579,6 +1579,83 @@ struct StrangenessInJetsIons { } } + /** + * @brief Add V0s as input for the jet finder algorithm in DATA + * + * @tparam Track + * @tparam V0PerColl The container type holding the sliced V0 candidates for the current collision. + * + * @param[in,out] fjInput Vector of FastJet PseudoJets where valid V0s will be appended. + * @param[in] fjTracks Vector containing tracks already selected for jet finder input. + * @param[in] v0sPerColl V0 candidates belonging to this specific collision. + * @param[in] vtxPos TVector3 object representing the vertex position. + */ + template + void AddV0sForJetReconstructionData(std::vector& fjInput, + const std::vector& fjTracks, + V0PerColl const& v0sPerColl, + const TVector3& vtxPos) + { + // Vector of labels: if true remove the element from fjInput + std::vector isTrackReplaced(fjTracks.size(), false); + std::vector v0PseudoJets; // V0s to use as input for jet finder + + for (const auto& v0 : v0sPerColl) { + const auto& pos = v0.template posTrack_as(); + const auto& neg = v0.template negTrack_as(); + + bool isK0S = false, isLambda = false, isAntiLambda = false; + // K0s + if (passedK0ShortSelection(v0, pos, neg, vtxPos)) { + // LOG(info) << "[AddV0sForJetReconstructionData] Add K0S as input for jet finder."; + double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassK0Short); + fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); + v0PseudoJets.emplace_back(fourMomentum); + isK0S = true; + } + // Lambda + if (passedLambdaSelection(v0, pos, neg, vtxPos)) { + // LOG(info) << "[AddV0sForJetReconstructionData] Add Lambda as input for jet finder."; + double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassLambda0); + fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); + v0PseudoJets.emplace_back(fourMomentum); + isLambda = true; + } + // AntiLambda + if (passedAntiLambdaSelection(v0, pos, neg, vtxPos)) { + // LOG(info) << "[AddV0sForJetReconstructionData] Add AntiLambda as input for jet finder."; + double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassLambda0Bar); + fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); + v0PseudoJets.emplace_back(fourMomentum); + isAntiLambda = true; + } + + // Exclude daughter tracks in case a V0 is found + bool isV0 = isK0S || isLambda || isAntiLambda; + if (!isV0) + continue; + for (int i = 0; i < fjTracks.size(); ++i) { + if (isV0DaughterTrack(fjTracks[i], v0)) { + // LOG(info) << "[AddV0sForJetReconstructionData] V0 daughter track found in fjTracks."; + isTrackReplaced[i] = true; + } + } + } + + std::vector cleanFjInput; + cleanFjInput.reserve(fjInput.size()); + for (size_t i = 0; i < fjInput.size(); ++i) { + if (!isTrackReplaced[i]) + cleanFjInput.push_back(fjInput[i]); + } + for (const auto& v0pj : v0PseudoJets) { + cleanFjInput.push_back(v0pj); + } + // LOG(info) << "[AddV0sForJetReconstructionData] Size fjInput: " << fjInput.size(); + fjInput = std::move(cleanFjInput); + // LOG(info) << "[AddV0sForJetReconstructionData] Size fjInput dopo move: " << fjInput.size(); + } + /** * @brief Add V0s as input for the jet finder algorithm in MC reco (detector level) * @@ -1592,11 +1669,11 @@ struct StrangenessInJetsIons { * @param[in] vtxPos TVector3 object representing the vertex position. */ template - void AddV0sForJetReconstructionMC(std::vector& fjInput, - const std::vector& fjTracks, - V0PerColl const& v0sPerColl, - aod::McParticles const& mcParticles, - const TVector3& vtxPos) + void AddV0sForJetReconstructionMCD(std::vector& fjInput, + const std::vector& fjTracks, + V0PerColl const& v0sPerColl, + aod::McParticles const& mcParticles, + const TVector3& vtxPos) { // Vector of labels: if true remove the element from fjInput std::vector isTrackReplaced(fjTracks.size(), false); @@ -1628,7 +1705,7 @@ struct StrangenessInJetsIons { bool isK0S = false, isLambda = false, isAntiLambda = false; // K0s if (passedK0ShortSelection(v0, pos, neg, vtxPos) && pdgParent == kK0Short) { - // LOG(info) << "[AddV0sForJetReconstructionMC] Add K0S as input for jet finder."; + // LOG(info) << "[AddV0sForJetReconstructionMCD] Add K0S as input for jet finder."; double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassK0Short); fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); v0PseudoJets.emplace_back(fourMomentum); @@ -1636,7 +1713,7 @@ struct StrangenessInJetsIons { } // Lambda if (passedLambdaSelection(v0, pos, neg, vtxPos) && pdgParent == kLambda0) { - // LOG(info) << "[AddV0sForJetReconstructionMC] Add Lambda as input for jet finder."; + // LOG(info) << "[AddV0sForJetReconstructionMCD] Add Lambda as input for jet finder."; double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassLambda0); fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); v0PseudoJets.emplace_back(fourMomentum); @@ -1644,7 +1721,7 @@ struct StrangenessInJetsIons { } // AntiLambda if (passedAntiLambdaSelection(v0, pos, neg, vtxPos) && pdgParent == kLambda0Bar) { - // LOG(info) << "[AddV0sForJetReconstructionMC] Add AntiLambda as input for jet finder."; + // LOG(info) << "[AddV0sForJetReconstructionMCD] Add AntiLambda as input for jet finder."; double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassLambda0Bar); fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); v0PseudoJets.emplace_back(fourMomentum); @@ -1657,7 +1734,7 @@ struct StrangenessInJetsIons { continue; for (int i = 0; i < fjTracks.size(); ++i) { if (isV0DaughterTrack(fjTracks[i], v0)) { - // LOG(info) << "[AddV0sForJetReconstructionMC] V0 daughter track found in fjTracks."; + // LOG(info) << "[AddV0sForJetReconstructionMCD] V0 daughter track found in fjTracks."; isTrackReplaced[i] = true; } } @@ -1672,9 +1749,9 @@ struct StrangenessInJetsIons { for (const auto& v0pj : v0PseudoJets) { cleanFjInput.push_back(v0pj); } - // LOG(info) << "[AddV0sForJetReconstructionMC] Size fjInput: " << fjInput.size(); + // LOG(info) << "[AddV0sForJetReconstructionMCD] Size fjInput: " << fjInput.size(); fjInput = std::move(cleanFjInput); - // LOG(info) << "[AddV0sForJetReconstructionMC] Size fjInput dopo move: " << fjInput.size(); + // LOG(info) << "[AddV0sForJetReconstructionMCD] Size fjInput dopo move: " << fjInput.size(); } // Process data @@ -1682,6 +1759,9 @@ struct StrangenessInJetsIons { aod::CascDataExt const& Cascades, DaughterTracks const& tracks, aod::BCsWithTimestamps const&) { + // Vertex position vector + TVector3 vtxPos(collision.posX(), collision.posY(), collision.posZ()); + // Fill event counter before event selection registryData.fill(HIST("number_of_events_data"), 0.5); @@ -1738,6 +1818,7 @@ struct StrangenessInJetsIons { // Loop over reconstructed tracks std::vector fjParticles; + std::vector> fjTracks; for (auto const& track : tracks) { // Require that tracks pass selection criteria @@ -1747,7 +1828,14 @@ struct StrangenessInJetsIons { // 4-momentum representation of a particle fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(o2::constants::physics::MassPionCharged)); fjParticles.emplace_back(fourMomentum); + fjTracks.push_back(track); + } + + // Include V0s as tracks for jet reconstruction + if (useV0inJetRec && particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + AddV0sForJetReconstructionData(fjParticles, fjTracks, fullV0s, vtxPos); } + fjTracks.clear(); // Reject empty events if (fjParticles.size() < 1) @@ -1839,9 +1927,6 @@ struct StrangenessInJetsIons { const float deltaPhiUe2 = getDeltaPhi(v0dir.Phi(), ue2[i].Phi()); const float deltaRue2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2); - // Vertex position vector - TVector3 vtxPos(collision.posX(), collision.posY(), collision.posZ()); - // Fill Armenteros-Podolanski TH2 // registryQC.fill(HIST("ArmenterosPreSel_DATA"), v0.alpha(), v0.qtarm()); @@ -2406,7 +2491,7 @@ struct StrangenessInJetsIons { // Include V0s as tracks for jet reconstruction if (useV0inJetRec && particleOfInterestDict[ParticleOfInterest::kV0Particles]) { - AddV0sForJetReconstructionMC(fjParticles, fjTracks, v0sPerColl, mcParticles, vtxPos); + AddV0sForJetReconstructionMCD(fjParticles, fjTracks, v0sPerColl, mcParticles, vtxPos); } fjTracks.clear(); From e8a5192c28c4a2701c2f90c9eca3fcc116f06ee7 Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Mon, 18 May 2026 11:39:24 +0200 Subject: [PATCH 5/9] Add V0s in jet reconstruction in MC GEN with configurable useV0inJetRec --- .../Strangeness/strangenessInJetsIons.cxx | 106 +++++++++++++++--- 1 file changed, 88 insertions(+), 18 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index 72869b1f5e2..f2eba1e855b 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -1644,7 +1644,7 @@ struct StrangenessInJetsIons { std::vector cleanFjInput; cleanFjInput.reserve(fjInput.size()); - for (size_t i = 0; i < fjInput.size(); ++i) { + for (int i = 0; i < fjInput.size(); ++i) { if (!isTrackReplaced[i]) cleanFjInput.push_back(fjInput[i]); } @@ -1698,13 +1698,13 @@ struct StrangenessInJetsIons { if (motherId < 0) continue; // if motherId is -1, it means no common mother was found const auto& mother = mcParticles.iteratorAt(motherId); - // const bool isPhysPrim = mother.isPhysicalPrimary(); + const bool isPhysPrim = mother.isPhysicalPrimary(); int pdgParent = mother.pdgCode(); - /// NOTE: V0s are not required to be primary particles; + /// NOTE: V0s are required to be primary particles; bool isK0S = false, isLambda = false, isAntiLambda = false; // K0s - if (passedK0ShortSelection(v0, pos, neg, vtxPos) && pdgParent == kK0Short) { + if (passedK0ShortSelection(v0, pos, neg, vtxPos) && pdgParent == kK0Short && isPhysPrim) { // LOG(info) << "[AddV0sForJetReconstructionMCD] Add K0S as input for jet finder."; double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassK0Short); fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); @@ -1712,7 +1712,7 @@ struct StrangenessInJetsIons { isK0S = true; } // Lambda - if (passedLambdaSelection(v0, pos, neg, vtxPos) && pdgParent == kLambda0) { + if (passedLambdaSelection(v0, pos, neg, vtxPos) && pdgParent == kLambda0 && isPhysPrim) { // LOG(info) << "[AddV0sForJetReconstructionMCD] Add Lambda as input for jet finder."; double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassLambda0); fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); @@ -1720,7 +1720,7 @@ struct StrangenessInJetsIons { isLambda = true; } // AntiLambda - if (passedAntiLambdaSelection(v0, pos, neg, vtxPos) && pdgParent == kLambda0Bar) { + if (passedAntiLambdaSelection(v0, pos, neg, vtxPos) && pdgParent == kLambda0Bar && isPhysPrim) { // LOG(info) << "[AddV0sForJetReconstructionMCD] Add AntiLambda as input for jet finder."; double energy = GetEnergy(v0.px(), v0.py(), v0.pz(), o2::constants::physics::MassLambda0Bar); fastjet::PseudoJet fourMomentum(v0.px(), v0.py(), v0.pz(), energy); @@ -1742,7 +1742,7 @@ struct StrangenessInJetsIons { std::vector cleanFjInput; cleanFjInput.reserve(fjInput.size()); - for (size_t i = 0; i < fjInput.size(); ++i) { + for (int i = 0; i < fjInput.size(); ++i) { if (!isTrackReplaced[i]) cleanFjInput.push_back(fjInput[i]); } @@ -1754,6 +1754,68 @@ struct StrangenessInJetsIons { // LOG(info) << "[AddV0sForJetReconstructionMCD] Size fjInput dopo move: " << fjInput.size(); } + /** + * @brief Add V0s as input for the jet finder algorithm at Particle Level (MC Generated) + * + * @tparam McParticleType + * @tparam ParticleColl The container type holding the sliced MC particles. + * + * @param[in,out] fjInput Vector of FastJet PseudoJets where V0s will be appended. + * @param[in] fjParticleObj Vector containing the MC particles already selected for jet finder input. + * @param[in] mcParticlesPerColl MC particles belonging to this specific collision. + * @param[in] mcParticles Full McParticles. + */ + template + void AddV0sForJetReconstructionMCP(std::vector& fjInput, + const std::vector& fjParticleObj, + ParticleColl const& mcParticlesPerColl, + aod::McParticles const& mcParticles) + { + std::vector isTrackReplaced(fjParticleObj.size(), false); + std::vector v0PseudoJets; + + // Add V0s to the input list for the jet finder and eventually remove their daughter particles + for (const auto& p : mcParticlesPerColl) { + int pdgAbs = std::abs(p.pdgCode()); + /// NOTE: p.isPhysicalPrimary() is required + if (p.isPhysicalPrimary() && (pdgAbs == kK0Short || pdgAbs == kLambda0)) { + double mass = (pdgAbs == kK0Short) ? o2::constants::physics::MassK0Short : o2::constants::physics::MassLambda0; + double energy = std::sqrt(p.p() * p.p() + mass * mass); + + fastjet::PseudoJet pj(p.px(), p.py(), p.pz(), energy); + pj.set_user_index(p.pdgCode()); + v0PseudoJets.push_back(pj); + // LOG(info) << "[AddV0sForJetReconstructionMCP] Add V0 as input for jet finder."; + + // Remove V0 daughter particles if already in the input list for the jet finder + for (int i = 0; i < fjParticleObj.size(); ++i) { + const auto& mcPart = fjParticleObj[i]; + if (!mcPart.has_mothers()) + continue; + auto mother = mcParticles.iteratorAt(mcPart.mothersIds()[0]); + int motherPdg = std::abs(mother.pdgCode()); + if (motherPdg == kK0Short || motherPdg == kLambda0) { + isTrackReplaced[i] = true; + // LOG(info) << "[AddV0sForJetReconstructionMCP] V0 daughter particle found in fjParticleObj."; + } + } + } + } + + std::vector cleanFjInput; + cleanFjInput.reserve(fjInput.size()); + for (int i = 0; i < fjInput.size(); ++i) { + if (!isTrackReplaced[i]) + cleanFjInput.push_back(fjInput[i]); + } + for (const auto& v0pj : v0PseudoJets) { + cleanFjInput.push_back(v0pj); + } + // LOG(info) << "[AddV0sForJetReconstructionMCP] Size fjInput: " << fjInput.size(); + fjInput = std::move(cleanFjInput); + // LOG(info) << "[AddV0sForJetReconstructionMCP] Size fjInput dopo move: " << fjInput.size(); + } + // Process data void processData(SelCollisions::iterator const& collision, aod::V0Datas const& fullV0s, aod::CascDataExt const& Cascades, DaughterTracks const& tracks, @@ -2158,6 +2220,7 @@ struct StrangenessInJetsIons { auto mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, collision.globalIndex()); // Loop over all MC particles and select physical primaries within acceptance + std::vector> fjParticleObj; for (const auto& particle : mcParticlesPerColl) { // Store properties of strange hadrons int pdgAbs = std::abs(particle.pdgCode()); @@ -2182,6 +2245,11 @@ struct StrangenessInJetsIons { fastjet::PseudoJet fourMomentum(particle.px(), particle.py(), particle.pz(), energy); fourMomentum.set_user_index(particle.pdgCode()); fjParticles.emplace_back(fourMomentum); + fjParticleObj.push_back(particle); + } + + if (useV0inJetRec && particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + AddV0sForJetReconstructionMCP(fjParticles, fjParticleObj, mcParticlesPerColl, mcParticles); } // Skip events with no particles @@ -2650,14 +2718,14 @@ struct StrangenessInJetsIons { continue; // --- alternative method to avoid double for loop - int motherId = GetCommonMotherId(mcParticles, posParticle, negParticle); - if (motherId < 0) - continue; // if motherId is -1, it means no common mother was found - const auto& motherTest = mcParticles.iteratorAt(motherId); - // const bool isPhysPrimTest = motherTest.isPhysicalPrimary(); - int pdgParentTest = motherTest.pdgCode(); - if (pdgParent != pdgParentTest) - LOG(warning) << "pdgParent != pdgParentTest" << endl; + // int motherId = GetCommonMotherId(mcParticles, posParticle, negParticle); + // if (motherId < 0) + // continue; // if motherId is -1, it means no common mother was found + // const auto& motherTest = mcParticles.iteratorAt(motherId); + // // const bool isPhysPrimTest = motherTest.isPhysicalPrimary(); + // int pdgParentTest = motherTest.pdgCode(); + // if (pdgParent != pdgParentTest) + // LOG(warning) << "pdgParent != pdgParentTest" << endl; // --- // Compute distance from jet and UE axes @@ -2854,9 +2922,6 @@ struct StrangenessInJetsIons { // Fill event counter after selection on z-vertex registryMC.fill(HIST("number_of_events_mc_gen"), 1.5); - // Multiplicity of generated event retrived from corresponding MC RECO collision - // float genMultiplicity = multiplicity; - // Multiplicity of generated event float genMultiplicity; if (centrEstimator == 0) { @@ -2870,6 +2935,7 @@ struct StrangenessInJetsIons { auto mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, collision.globalIndex()); // Loop over all MC particles and select physical primaries within acceptance + std::vector> fjParticleObj; for (const auto& particle : mcParticlesPerColl) { // Store properties of strange hadrons int pdgAbs = std::abs(particle.pdgCode()); @@ -2897,6 +2963,10 @@ struct StrangenessInJetsIons { fjParticles.emplace_back(fourMomentum); } + if (useV0inJetRec && particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + AddV0sForJetReconstructionMCP(fjParticles, fjParticleObj, mcParticlesPerColl, mcParticles); + } + // Skip events with no particles if (fjParticles.size() < 1) continue; From 9081756e50073afcfea1668c56a19cb06685bbad Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Mon, 18 May 2026 11:49:15 +0200 Subject: [PATCH 6/9] Fix bug in processMCmodels: fill fjParticleObj --- PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index f2eba1e855b..7dac0b43df5 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -2961,6 +2961,7 @@ struct StrangenessInJetsIons { fastjet::PseudoJet fourMomentum(particle.px(), particle.py(), particle.pz(), energy); fourMomentum.set_user_index(particle.pdgCode()); fjParticles.emplace_back(fourMomentum); + fjParticleObj.push_back(particle); } if (useV0inJetRec && particleOfInterestDict[ParticleOfInterest::kV0Particles]) { From c2fdbaa7d87b312d0498260847acbc770af07877 Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Mon, 18 May 2026 12:15:31 +0200 Subject: [PATCH 7/9] Add include utility for std::move --- PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index 7dac0b43df5..ca795e53e8f 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -71,6 +71,7 @@ #include #include #include +#include using namespace std; using namespace o2; From 13e4a4bfc8a23a1c954c6e756ec11d90a38db975 Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Mon, 18 May 2026 12:25:20 +0200 Subject: [PATCH 8/9] Address formatting errors --- PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index ca795e53e8f..12dba4f3e24 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -22,7 +22,6 @@ #include "PWGJE/Core/JetBkgSubUtils.h" #include "PWGJE/Core/JetUtilities.h" -// #include "PWGJE/Core/JetV0Utilities.h" #include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGLF/DataModel/mcCentrality.h" @@ -70,8 +69,8 @@ #include #include #include -#include #include +#include using namespace std; using namespace o2; From a85c2e423e66dc2713b49c71f5bc86168dbba77c Mon Sep 17 00:00:00 2001 From: lorber98 <95907752+lorber98@users.noreply.github.com> Date: Mon, 18 May 2026 13:01:53 +0200 Subject: [PATCH 9/9] Fix comparison fjInput.size() in for loops --- PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index 12dba4f3e24..6d654706c2e 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -1634,7 +1634,7 @@ struct StrangenessInJetsIons { bool isV0 = isK0S || isLambda || isAntiLambda; if (!isV0) continue; - for (int i = 0; i < fjTracks.size(); ++i) { + for (int i = 0; i < int(fjTracks.size()); ++i) { if (isV0DaughterTrack(fjTracks[i], v0)) { // LOG(info) << "[AddV0sForJetReconstructionData] V0 daughter track found in fjTracks."; isTrackReplaced[i] = true; @@ -1644,7 +1644,7 @@ struct StrangenessInJetsIons { std::vector cleanFjInput; cleanFjInput.reserve(fjInput.size()); - for (int i = 0; i < fjInput.size(); ++i) { + for (int i = 0; i < int(fjInput.size()); ++i) { if (!isTrackReplaced[i]) cleanFjInput.push_back(fjInput[i]); } @@ -1732,7 +1732,7 @@ struct StrangenessInJetsIons { bool isV0 = isK0S || isLambda || isAntiLambda; if (!isV0) continue; - for (int i = 0; i < fjTracks.size(); ++i) { + for (int i = 0; i < int(fjTracks.size()); ++i) { if (isV0DaughterTrack(fjTracks[i], v0)) { // LOG(info) << "[AddV0sForJetReconstructionMCD] V0 daughter track found in fjTracks."; isTrackReplaced[i] = true; @@ -1742,7 +1742,7 @@ struct StrangenessInJetsIons { std::vector cleanFjInput; cleanFjInput.reserve(fjInput.size()); - for (int i = 0; i < fjInput.size(); ++i) { + for (int i = 0; i < int(fjInput.size()); ++i) { if (!isTrackReplaced[i]) cleanFjInput.push_back(fjInput[i]); } @@ -1788,7 +1788,7 @@ struct StrangenessInJetsIons { // LOG(info) << "[AddV0sForJetReconstructionMCP] Add V0 as input for jet finder."; // Remove V0 daughter particles if already in the input list for the jet finder - for (int i = 0; i < fjParticleObj.size(); ++i) { + for (int i = 0; i < int(fjParticleObj.size()); ++i) { const auto& mcPart = fjParticleObj[i]; if (!mcPart.has_mothers()) continue; @@ -1804,7 +1804,7 @@ struct StrangenessInJetsIons { std::vector cleanFjInput; cleanFjInput.reserve(fjInput.size()); - for (int i = 0; i < fjInput.size(); ++i) { + for (int i = 0; i < int(fjInput.size()); ++i) { if (!isTrackReplaced[i]) cleanFjInput.push_back(fjInput[i]); }