diff --git a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx index ffb3b001f4d..6d654706c2e 100644 --- a/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx +++ b/PWGLF/Tasks/Strangeness/strangenessInJetsIons.cxx @@ -68,6 +68,8 @@ #include #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"}; @@ -334,7 +337,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"}}); @@ -346,9 +352,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}); @@ -558,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) { @@ -1187,9 +1201,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 +1289,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]) { @@ -1524,11 +1538,292 @@ 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 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 < int(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 (int i = 0; i < int(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) + * + * @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 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); + 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 required to be primary particles; + bool isK0S = false, isLambda = false, isAntiLambda = false; + // K0s + 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); + v0PseudoJets.emplace_back(fourMomentum); + isK0S = true; + } + // Lambda + 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); + v0PseudoJets.emplace_back(fourMomentum); + isLambda = true; + } + // AntiLambda + 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); + 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 < int(fjTracks.size()); ++i) { + if (isV0DaughterTrack(fjTracks[i], v0)) { + // LOG(info) << "[AddV0sForJetReconstructionMCD] V0 daughter track found in fjTracks."; + isTrackReplaced[i] = true; + } + } + } + + std::vector cleanFjInput; + cleanFjInput.reserve(fjInput.size()); + for (int i = 0; i < int(fjInput.size()); ++i) { + if (!isTrackReplaced[i]) + cleanFjInput.push_back(fjInput[i]); + } + for (const auto& v0pj : v0PseudoJets) { + cleanFjInput.push_back(v0pj); + } + // LOG(info) << "[AddV0sForJetReconstructionMCD] Size fjInput: " << fjInput.size(); + fjInput = std::move(cleanFjInput); + // 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 < int(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 < int(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, 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); @@ -1585,6 +1880,7 @@ struct StrangenessInJetsIons { // Loop over reconstructed tracks std::vector fjParticles; + std::vector> fjTracks; for (auto const& track : tracks) { // Require that tracks pass selection criteria @@ -1594,7 +1890,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) @@ -1686,9 +1989,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()); @@ -1913,11 +2213,14 @@ 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()); // 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()); @@ -1930,7 +2233,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)) @@ -1942,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 @@ -1974,6 +2282,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()); @@ -2181,6 +2491,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(); @@ -2229,10 +2542,11 @@ struct StrangenessInJetsIons { auto tracksPerColl = mcTracks.sliceBy(perCollisionTrk, collision.globalIndex()); const auto& mcParticlesPerColl = mcParticles.sliceBy(perMCCollision, mcCollision.globalIndex()); - FillFullEventHistoMCREC(collision, mcParticles, v0sPerColl, - cascPerColl, tracksPerColl, multiplicity); + FillMBEventHistoMCREC(collision, mcParticles, v0sPerColl, + cascPerColl, tracksPerColl, multiplicity); // Loop over reconstructed tracks + std::vector> fjTracks; for (auto const& track : tracksPerColl) { if (!passedTrackSelectionForJetReconstruction(track)) continue; @@ -2240,7 +2554,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]) { + AddV0sForJetReconstructionMCD(fjParticles, fjTracks, v0sPerColl, mcParticles, vtxPos); } + fjTracks.clear(); // Reject empty events if (fjParticles.size() < 1) @@ -2396,6 +2717,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()); @@ -2407,9 +2739,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()); @@ -2563,6 +2892,294 @@ 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 + 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 + std::vector> fjParticleObj; + 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); + fjParticleObj.push_back(particle); + } + + if (useV0inJetRec && particleOfInterestDict[ParticleOfInterest::kV0Particles]) { + AddV0sForJetReconstructionMCP(fjParticles, fjParticleObj, mcParticlesPerColl, mcParticles); + } + + // 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)