diff --git a/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx b/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx index 59aca0f7ceb..ad51876259b 100644 --- a/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx +++ b/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx @@ -77,6 +77,7 @@ using std::array; using SelectedCollisions = soa::Join; using SimCollisions = soa::Join; +using GenCollisions = aod::McCollisions; using FullNucleiTracks = soa::Join; @@ -88,9 +89,11 @@ struct AntinucleiInJets { HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry registryMC{"registryMC", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry registryQC{"registryQC", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry registryMult{"registryMult", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; // global parameters Configurable minJetPt{"minJetPt", 10.0, "Minimum pt of the jet"}; + Configurable ptLeadingMin{"ptLeadingMin", 5.0, "pt Leading Min"}; Configurable rJet{"rJet", 0.3, "Jet resolution parameter R"}; Configurable zVtx{"zVtx", 10.0, "Maximum zVertex"}; Configurable deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"}; @@ -179,6 +182,11 @@ struct AntinucleiInJets { registryQC.add("jetPtDifference", "jetPtDifference", HistType::kTH1F, {{200, -1, 1, "#Deltap_{T}^{jet}"}}); } + if (doprocessMultEvents) { + registryMult.add("multiplicityEvtsPtLeading", "multiplicityEvtsPtLeading", HistType::kTH1F, {{1000, 0, 1000, "#it{N}_{ch}"}}); + registryMult.add("multiplicityEvtsWithJet", "multiplicityEvtsWithJet", HistType::kTH1F, {{1000, 0, 1000, "#it{N}_{ch}"}}); + } + // data histograms if (doprocessData) { @@ -191,8 +199,8 @@ struct AntinucleiInJets { registryData.add("antiproton_jet_tof", "antiproton_jet_tof", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TOF}"}}); registryData.add("antiproton_ue_tpc", "antiproton_ue_tpc", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TPC}"}}); registryData.add("antiproton_ue_tof", "antiproton_ue_tof", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TOF}"}}); - registryData.add("antiproton_dca_jet", "antiproton_dca_jet", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -0.5, 0.5, "DCA_{xy} (cm)"}}); - registryData.add("antiproton_dca_ue", "antiproton_dca_ue", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -0.5, 0.5, "DCA_{xy} (cm)"}}); + registryData.add("antiproton_dca_jet", "antiproton_dca_jet", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); + registryData.add("antiproton_dca_ue", "antiproton_dca_ue", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); // antideuterons registryData.add("antideuteron_jet_tpc", "antideuteron_jet_tpc", HistType::kTH2F, {{nbins, min * 2, max * 2, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{TPC}"}}); @@ -271,6 +279,12 @@ struct AntinucleiInJets { registryMC.add("antiproton_prim_ue", "antiproton_prim_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); registryMC.add("antiproton_incl_ue", "antiproton_incl_ue", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); + // DCA Templates + registryData.add("antiproton_prim_dca_jet", "antiproton_prim_dca_jet", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); + registryData.add("antiproton_prim_dca_ue", "antiproton_prim_dca_ue", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); + registryData.add("antiproton_all_dca_jet", "antiproton_all_dca_jet", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); + registryData.add("antiproton_all_dca_ue", "antiproton_all_dca_ue", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {200, -1, 1, "DCA_{xy} (cm)"}}); + // antiproton reweighting registryMC.add("antiproton_eta_pt_pythia", "antiproton_eta_pt_pythia", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}}); registryMC.add("antiproton_eta_pt_jet", "antiproton_eta_pt_jet", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}}); @@ -597,7 +611,8 @@ struct AntinucleiInJets { // jet pt must be larger than threshold auto jetForSub = jet; fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp); - if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt) + // if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt) + if (jetMinusBkg.pt() < minJetPt) continue; isAtLeastOneJetSelected = true; @@ -768,6 +783,68 @@ struct AntinucleiInJets { } PROCESS_SWITCH(AntinucleiInJets, processData, "Process Data", true); + void processMultEvents(SelectedCollisions::iterator const& collision, FullNucleiTracks const& tracks) + { + // event selection + if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) + return; + + // Leading Track + double ptMax(0.0); + + // loop over reconstructed tracks + int id(-1); + std::vector fjParticles; + for (auto const& track : tracks) { + id++; + if (!passedTrackSelectionForJetReconstruction(track)) + continue; + if (track.pt() > ptMax) { + ptMax = track.pt(); + } + + // 4-momentum representation of a particle + fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged)); + fourMomentum.set_user_index(id); + fjParticles.emplace_back(fourMomentum); + } + // reject empty events + if (fjParticles.empty()) { + return; + } + + if (ptMax > ptLeadingMin) { + registryMult.fill(HIST("multiplicityEvtsPtLeading"), fjParticles.size()); + } + + // cluster particles using the anti-kt algorithm + fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet); + fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0)); + fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); + std::vector jets = fastjet::sorted_by_pt(cs.inclusive_jets()); + auto [rhoPerp, rhoMPerp] = backgroundSub.estimateRhoPerpCone(fjParticles, jets); + + // loop over reconstructed jets + bool isAtLeastOneJetSelected = false; + for (const auto& jet : jets) { + + // jet must be fully contained in the acceptance + if ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge)) + continue; + + // jet pt must be larger than threshold + auto jetForSub = jet; + fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp); + if (jetMinusBkg.pt() < minJetPt) + continue; + isAtLeastOneJetSelected = true; + } + if (isAtLeastOneJetSelected) { + registryMult.fill(HIST("multiplicityEvtsWithJet"), fjParticles.size()); + } + } + PROCESS_SWITCH(AntinucleiInJets, processMultEvents, "Process Mult Events", true); + // Process QC void processQC(SelectedCollisions::iterator const& collision, FullNucleiTracks const& tracks) { @@ -1087,13 +1164,13 @@ struct AntinucleiInJets { } PROCESS_SWITCH(AntinucleiInJets, processEfficiency, "process efficiency", false); - void processJetsMCgen(SimCollisions const& collisions, aod::McParticles const& mcParticles) + void processJetsMCgen(GenCollisions const& collisions, aod::McParticles const& mcParticles) { // Loop over all simulated collision events for (const auto& collision : collisions) { - // Apply event selection: require sel8 and vertex position within the allowed z range - if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) + // Apply event selection: require vertex position within the allowed z range + if (std::fabs(collision.posZ()) > zVtx) continue; // Loop over all MC particles and select physical primaries within acceptance @@ -1256,7 +1333,8 @@ struct AntinucleiInJets { fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp); // Apply jet pT threshold - if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt) + // if (getCorrectedPt(jetMinusBkg.pt(), responseMatrix) < minJetPt) + if (jetMinusBkg.pt() < minJetPt) continue; // Set up two perpendicular cone axes for underlying event estimation @@ -1271,11 +1349,20 @@ struct AntinucleiInJets { auto const& track = mcTracks.iteratorAt(particle.user_index()); if (!track.has_mcParticle()) continue; + const auto mcparticle = track.mcParticle(); + + // Fill DCA templates + if (mcparticle.pdgCode() == kProtonBar && passedTrackSelection(track) && std::fabs(track.dcaZ()) < maxDcaz) { + if (mcparticle.isPhysicalPrimary()) { + registryMC.fill(HIST("antiproton_prim_dca_jet"), track.pt(), track.dcaXY()); + } else { + registryMC.fill(HIST("antiproton_all_dca_jet"), track.pt(), track.dcaXY()); + } + } // Apply standard track quality and PID selection if (!passedTrackSelection(track) || std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; - const auto mcparticle = track.mcParticle(); if (track.sign() > 0 || mcparticle.pdgCode() != kProtonBar) continue; @@ -1310,15 +1397,26 @@ struct AntinucleiInJets { // Analyze antiprotons in the Underlying Event (UE) using perpendicular cones for (auto const& track : mcTracks) { - if (!passedTrackSelection(track) || std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) - continue; + if (!track.has_mcParticle()) continue; - const auto mcparticle = track.mcParticle(); if (track.sign() > 0 || mcparticle.pdgCode() != kProtonBar) continue; + // Fill DCA templates + if (mcparticle.pdgCode() == kProtonBar && passedTrackSelection(track) && std::fabs(track.dcaZ()) < maxDcaz) { + if (mcparticle.isPhysicalPrimary()) { + registryMC.fill(HIST("antiproton_prim_dca_ue"), track.pt(), track.dcaXY()); + } else { + registryMC.fill(HIST("antiproton_all_dca_ue"), track.pt(), track.dcaXY()); + } + } + + // Track Selection + if (!passedTrackSelection(track) || std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) + continue; + // Compute distance from UE cones double deltaEtaUe1 = track.eta() - ueAxis1.Eta(); double deltaPhiUe1 = getDeltaPhi(track.phi(), ueAxis1.Phi());