From f510ec1fbb3a13ef25f06377ff8edd126d90c7b2 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Mon, 12 Oct 2026 18:17:07 +0900 Subject: [PATCH 1/2] [PWGJE] Add baryon pair correlations in charged jets --- PWGJE/Tasks/CMakeLists.txt | 4 + PWGJE/Tasks/baryonPairInJets.cxx | 1110 ++++++++++++++++++++++++++++++ 2 files changed, 1114 insertions(+) create mode 100644 PWGJE/Tasks/baryonPairInJets.cxx diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index c633110950f..e010858ccbe 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -519,4 +519,8 @@ if(FastJet_FOUND) SOURCES embeddingQA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(baryon-pair-in-jets + SOURCES baryonPairInJets.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore + COMPONENT_NAME Analysis) endif() diff --git a/PWGJE/Tasks/baryonPairInJets.cxx b/PWGJE/Tasks/baryonPairInJets.cxx new file mode 100644 index 00000000000..bddcf690ae7 --- /dev/null +++ b/PWGJE/Tasks/baryonPairInJets.cxx @@ -0,0 +1,1110 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file baryonPairInJets.cxx +/// \brief Identified proton-pair correlations in charged jets +/// \author Hanseo Park + +#include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/DataModel/Jet.h" +#include "PWGJE/DataModel/JetReducedData.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; + +using BaryonTracks = soa::Join; +using BaryonTracksMC = soa::Join; +using JetTracksWithOriginal = soa::Join; +using JetTracksMCDWithOriginal = soa::Join; +using ChargedJetsWithConstituents = soa::Join; +using ChargedMCDJetsWithConstituents = soa::Join; +using ChargedMCPJetsWithConstituents = soa::Join; +using MatchedChargedMCDJets = soa::Join; +using MatchedChargedMCPJets = soa::Join; + +struct BaryonPairInJets { + enum class PairType : uint8_t { + SameSign = 0, + OppositeSign = 1 + }; + + enum class ChargePairType : uint8_t { + ProtonProton = 0, + AntiprotonAntiproton = 1, + ProtonAntiproton = 2 + }; + + enum class MixedReference : uint8_t { + Nominal = 0, + Loose = 1, + Random = 2 + }; + + enum class ProtonPidMethod : int { + TpcOnly = 0, + TpcTofCircular = 1, + TpcTofRectangular = 2 + }; + + struct SelectedBaryon { + float pt; + float eta; + float phi; + float rapidity; + int8_t sign; + }; + + struct : ConfigurableGroup { + Configurable vertexZMax{"vertexZMax", 10.f, "Maximum absolute collision z vertex (cm)"}; + Configurable eventSelections{"eventSelections", "sel8+NoITSROFrameBorder+NoTimeFrameBorder+NoSameBunchPileup+IsGoodZvtxFT0vsPV", "Event selections"}; + Configurable applyRctSelection{"applyRctSelection", true, "Apply the PWGJE hadron-PID RCT selection"}; + Configurable occupancyMin{"occupancyMin", -999999, "Minimum track occupancy in time range"}; + Configurable occupancyMax{"occupancyMax", 999999, "Maximum track occupancy in time range"}; + } eventCuts; + + struct : ConfigurableGroup { + Configurable jetRadius{"jetRadius", 0.4f, "Charged-jet resolution parameter"}; + Configurable jetPtMin{"jetPtMin", 5.f, "Minimum charged-jet transverse momentum (GeV/c); broad input range"}; + Configurable jetPtMax{"jetPtMax", 140.f, "Maximum charged-jet transverse momentum (GeV/c); broad input range"}; + Configurable jetEtaMax{"jetEtaMax", 0.5f, "Maximum absolute charged-jet pseudorapidity"}; + Configurable areaFractionMin{"areaFractionMin", -1.f, "Minimum jet area divided by pi R squared; negative disables the cut"}; + } jetCuts; + + struct : ConfigurableGroup { + Configurable protonPtMin{"protonPtMin", 0.5f, "Minimum proton transverse momentum (GeV/c)"}; + Configurable protonPtMax{"protonPtMax", 2.5f, "Maximum proton transverse momentum (GeV/c)"}; + Configurable protonRapidityMax{"protonRapidityMax", 0.5f, "Maximum absolute proton rapidity"}; + Configurable protonEtaMax{"protonEtaMax", 0.8f, "Maximum absolute proton pseudorapidity before the rapidity cut"}; + Configurable requirePvContributor{"requirePvContributor", false, "Require the track to be a PV contributor"}; + Configurable tpcCrossedRowsMin{"tpcCrossedRowsMin", 70, "Minimum number of crossed TPC rows"}; + Configurable tpcChi2Max{"tpcChi2Max", 4.f, "Maximum TPC chi2 per cluster"}; + Configurable itsChi2Max{"itsChi2Max", 36.f, "Maximum ITS chi2 per cluster"}; + Configurable dcaXYMax{"dcaXYMax", 0.2f, "Maximum absolute transverse DCA (cm)"}; + Configurable dcaZMax{"dcaZMax", 0.1f, "Maximum absolute longitudinal DCA (cm)"}; + } trackCuts; + + struct : ConfigurableGroup { + Configurable method{"method", static_cast(ProtonPidMethod::TpcTofCircular), "PID method: 0 TPC only, 1 circular TPC-TOF, 2 rectangular TPC-TOF"}; + Configurable tofPtThreshold{"tofPtThreshold", 0.75f, "Require TOF above this transverse momentum (GeV/c)"}; + Configurable tpcNSigmaMax{"tpcNSigmaMax", 3.f, "Maximum absolute proton TPC n-sigma"}; + Configurable tofNSigmaMax{"tofNSigmaMax", 3.f, "Maximum absolute proton TOF n-sigma"}; + Configurable combinedNSigmaMax{"combinedNSigmaMax", 3.f, "Maximum circular TPC-TOF proton n-sigma"}; + Configurable rejectPionKaon{"rejectPionKaon", false, "Reject candidates compatible with the pion or kaon hypothesis"}; + Configurable rejectionNSigma{"rejectionNSigma", 3.f, "Pion and kaon rejection threshold"}; + } pidCuts; + + struct : ConfigurableGroup { + Configurable numberOfMixedEvents{"numberOfMixedEvents", 5, "Number of mixed-event neighbours"}; + ConfigurableAxis vertexZBins{"vertexZBins", {VARIABLE_WIDTH, -10.f, -5.f, -2.f, 0.f, 2.f, 5.f, 10.f}, "Mixing bins in collision z vertex"}; + ConfigurableAxis multiplicityBins{"multiplicityBins", {VARIABLE_WIDTH, 0.f, 10.f, 20.f, 40.f, 80.f, 150.f, 300.f, 600.f, 2000.f}, "Mixing bins in FT0C multiplicity"}; + Configurable nominalJetPtDifferenceMax{"nominalJetPtDifferenceMax", 2.5f, "Maximum jet-pT difference for nominal mixing (GeV/c)"}; + Configurable nominalJetNchDifferenceMax{"nominalJetNchDifferenceMax", 3, "Maximum jet-constituent multiplicity difference for nominal mixing"}; + Configurable nominalJetEtaDifferenceMax{"nominalJetEtaDifferenceMax", 0.1f, "Maximum jet-eta difference for nominal mixing"}; + Configurable looseJetPtDifferenceMax{"looseJetPtDifferenceMax", 7.5f, "Maximum jet-pT difference for loose mixing (GeV/c)"}; + Configurable looseJetNchDifferenceMax{"looseJetNchDifferenceMax", 8, "Maximum jet-constituent multiplicity difference for loose mixing"}; + Configurable looseJetEtaDifferenceMax{"looseJetEtaDifferenceMax", 0.3f, "Maximum jet-eta difference for loose mixing"}; + } mixing; + + struct : ConfigurableGroup { + Configurable nearSideDeltaPhiMax{"nearSideDeltaPhiMax", 0.2f, "Maximum absolute delta phi for near-side pair yields"}; + Configurable nearSideDeltaEtaMax{"nearSideDeltaEtaMax", 0.6f, "Maximum absolute delta eta for near-side pair yields"}; + ConfigurableAxis jetPt{"jetPt", {VARIABLE_WIDTH, 0.f, 5.f, 10.f, 15.f, 20.f, 25.f, 30.f, 40.f, 50.f, 60.f, 80.f, 100.f, 120.f, 140.f}, "#it{p}_{T, jet} (GeV/#it{c})"}; + ConfigurableAxis jetEta{"jetEta", {40, -1.f, 1.f}, "#eta_{jet}"}; + ConfigurableAxis jetPhi{"jetPhi", {72, 0.f, o2::constants::math::TwoPI}, "#varphi_{jet}"}; + ConfigurableAxis jetNch{"jetNch", {61, -0.5f, 60.5f}, "#it{N}_{ch}^{jet}"}; + ConfigurableAxis baryonPt{"baryonPt", {50, 0.f, 5.f}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis baryonEta{"baryonEta", {48, -1.2f, 1.2f}, "#eta"}; + ConfigurableAxis baryonRapidity{"baryonRapidity", {48, -1.2f, 1.2f}, "#it{y}"}; + ConfigurableAxis nSigma{"nSigma", {100, -5.f, 5.f}, "n#sigma"}; + ConfigurableAxis dca{"dca", {100, -0.5f, 0.5f}, "DCA (cm)"}; + ConfigurableAxis deltaEta{"deltaEta", {80, -0.8f, 0.8f}, "#Delta#eta"}; + ConfigurableAxis deltaPhi{"deltaPhi", {80, -0.8f, 0.8f}, "#Delta#varphi"}; + ConfigurableAxis deltaR{"deltaR", {40, 0.f, 0.8f}, "#Delta#it{R}"}; + ConfigurableAxis inclusiveDeltaEta{"inclusiveDeltaEta", {64, -1.6f, 1.6f}, "#Delta#eta"}; + ConfigurableAxis inclusiveDeltaPhi{"inclusiveDeltaPhi", {72, -o2::constants::math::PIHalf, 3.f * o2::constants::math::PIHalf}, "#Delta#varphi"}; + } axes; + + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; + + std::vector eventSelectionBits; + SliceCache mixingCache; + SliceCache mcdMixingCache; + SliceCache mcpMixingCache; + SliceCache inclusiveMixingCache; + using MixingBinning = ColumnBinningPolicy; + using MCPMixingBinning = ColumnBinningPolicy; + using InclusiveCollisions = soa::Join; + using InclusiveMixingBinning = ColumnBinningPolicy; + + void init(InitContext const&) + { + eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventCuts.eventSelections)); + + AxisSpec axisPairType{2, -0.5, 1.5, "pair type"}; + AxisSpec axisChargePairType{3, -0.5, 2.5, "charge pair type"}; + AxisSpec axisCounter{8, 0.5, 8.5, "counter"}; + AxisSpec axisMultiplicity{11, -0.5, 10.5, "multiplicity per jet"}; + AxisSpec axisJetPt{axes.jetPt}; + AxisSpec axisJetEta{axes.jetEta}; + AxisSpec axisJetPhi{axes.jetPhi}; + AxisSpec axisJetNch{axes.jetNch}; + AxisSpec axisBaryonPt{axes.baryonPt}; + AxisSpec axisBaryonEta{axes.baryonEta}; + AxisSpec axisBaryonRapidity{axes.baryonRapidity}; + AxisSpec axisNSigma{axes.nSigma}; + AxisSpec axisDca{axes.dca}; + AxisSpec axisDeltaEta{axes.deltaEta}; + AxisSpec axisDeltaPhi{axes.deltaPhi}; + AxisSpec axisDeltaR{axes.deltaR}; + AxisSpec axisInclusiveDeltaEta{axes.inclusiveDeltaEta}; + AxisSpec axisInclusiveDeltaPhi{axes.inclusiveDeltaPhi}; + AxisSpec axisOrigin{3, 0.5, 3.5, "origin"}; + + registry.add("event/cutFlow", "Event selection;selection step;entries", HistType::kTH1F, {axisCounter}); + registry.add("event/zVertex", "Selected collision z vertex;#it{z}_{vtx} (cm);entries", HistType::kTH1F, {{80, -20.f, 20.f}}); + registry.add("event/multiplicityFT0C", "Selected collision FT0C multiplicity;FT0C multiplicity;entries", HistType::kTH1F, {{200, 0.f, 2000.f}}); + + registry.add("jet/ptBeforePtSelection", "Charged jets before the jet-pT selection;#it{p}_{T,jet} (GeV/#it{c});entries", HistType::kTH1F, {axisJetPt}); + registry.add("jet/pt", "Selected charged jets;#it{p}_{T,jet} (GeV/#it{c});entries", HistType::kTH1F, {axisJetPt}); + registry.add("jet/eta", "Selected charged jets;#eta_{jet};entries", HistType::kTH1F, {axisJetEta}); + registry.add("jet/phi", "Selected charged jets;#varphi_{jet};entries", HistType::kTH1F, {axisJetPhi}); + registry.add("jet/nConstituents", "Charged-jet constituent multiplicity;#it{N}_{ch}^{jet};entries", HistType::kTH1F, {axisJetNch}); + registry.add("jet/nSelectedBaryons", "Selected baryon multiplicity per jet;#it{N}_{p+#bar{p}};jets", HistType::kTH1F, {axisMultiplicity}); + registry.add("jet/nProtons", "Selected proton multiplicity per jet;#it{N}_{p};jets", HistType::kTH1F, {axisMultiplicity}); + registry.add("jet/nAntiprotons", "Selected antiproton multiplicity per jet;#it{N}_{#bar{p}};jets", HistType::kTH1F, {axisMultiplicity}); + registry.add("jet/nProtonsVsAntiprotons", "Proton and antiproton multiplicities per jet;#it{N}_{p};#it{N}_{#bar{p}}", HistType::kTH2F, {axisMultiplicity, axisMultiplicity}); + registry.add("jet/nConstituentsVsPt", "Constituent multiplicity versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});#it{N}_{ch}^{jet}", HistType::kTH2F, {axisJetPt, axisJetNch}); + registry.add("jet/nBaryonsVsPt", "Selected baryon multiplicity versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});#it{N}_{p+#bar{p}}", HistType::kTH2F, {axisJetPt, axisMultiplicity}); + registry.add("jet/sameSignPairsVsPt", "Same-sign pair count versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});pairs", HistType::kTH1F, {axisJetPt}); + registry.add("jet/oppositeSignPairsVsPt", "Opposite-sign pair count versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});pairs", HistType::kTH1F, {axisJetPt}); + + registry.add("proton/pt", "Selected proton and antiproton candidates;#it{p}_{T} (GeV/#it{c});entries", HistType::kTH1F, {axisBaryonPt}); + registry.add("proton/eta", "Selected proton and antiproton candidates;#eta;entries", HistType::kTH1F, {axisBaryonEta}); + registry.add("proton/rapidity", "Selected proton and antiproton candidates;#it{y};entries", HistType::kTH1F, {axisBaryonRapidity}); + registry.add("proton/tpcNSigmaVsPt", "Proton TPC PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("proton/tofNSigmaVsPt", "Proton TOF PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("proton/dcaXYVsPt", "Selected proton transverse DCA;#it{p}_{T} (GeV/#it{c});DCA_{xy} (cm)", HistType::kTH2F, {axisBaryonPt, axisDca}); + registry.add("proton/dcaZVsPt", "Selected proton longitudinal DCA;#it{p}_{T} (GeV/#it{c});DCA_{z} (cm)", HistType::kTH2F, {axisBaryonPt, axisDca}); + registry.add("pid/before/tpcSignalVsSignedP", "TPC response before proton PID;charge #times #it{p} (GeV/#it{c});TPC signal", HistType::kTH2F, {{120, -6.f, 6.f}, {200, 0.f, 200.f}}); + registry.add("pid/before/tofBetaVsSignedP", "TOF response before proton PID;charge #times #it{p} (GeV/#it{c});#beta", HistType::kTH2F, {{120, -6.f, 6.f}, {120, 0.f, 1.2f}}); + registry.add("pid/before/tpcNSigmaPiVsPt", "Pion TPC PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{#pi}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tpcNSigmaKaVsPt", "Kaon TPC PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{K}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tpcNSigmaPrVsPt", "Proton TPC PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tofNSigmaPiVsPt", "Pion TOF PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{#pi}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tofNSigmaKaVsPt", "Kaon TOF PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{K}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tofNSigmaPrVsPt", "Proton TOF PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/tofMatchDenominator", "Tracks entering PID QA;#it{p}_{T} (GeV/#it{c});tracks", HistType::kTH1F, {axisBaryonPt}); + registry.add("pid/tofMatchNumerator", "Tracks with TOF;#it{p}_{T} (GeV/#it{c});tracks", HistType::kTH1F, {axisBaryonPt}); + registry.add("pid/selectedPtByCharge", "Selected proton candidates;charge;#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{2, -1.5f, 1.5f}, axisBaryonPt}); + + addPairHistograms("sameEvent", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + addPairHistograms("mixedEvent/nominal", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + addPairHistograms("mixedEvent/loose", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + addPairHistograms("mixedEvent/random", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + registry.add("sameEvent/chargePairSparse", "Charge-separated baryon pairs;charge pair type;#it{p}_{T,jet};#Delta#eta;#Delta#varphi;#Delta#it{R}", HistType::kTHnSparseF, {axisChargePairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR}); + registry.add("mixedEvent/nominal/chargePairSparse", "Charge-separated nominal mixed pairs;charge pair type;#it{p}_{T,jet};#Delta#eta;#Delta#varphi;#Delta#it{R}", HistType::kTHnSparseF, {axisChargePairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR}); + registry.add("mixedEvent/loose/chargePairSparse", "Charge-separated loose mixed pairs;charge pair type;#it{p}_{T,jet};#Delta#eta;#Delta#varphi;#Delta#it{R}", HistType::kTHnSparseF, {axisChargePairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR}); + registry.add("mixedEvent/random/chargePairSparse", "Charge-separated random mixed pairs;charge pair type;#it{p}_{T,jet};#Delta#eta;#Delta#varphi;#Delta#it{R}", HistType::kTHnSparseF, {axisChargePairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR}); + + addInclusivePairHistograms("inclusive/sameEvent", axisChargePairType, axisInclusiveDeltaEta, axisInclusiveDeltaPhi); + addInclusivePairHistograms("inclusive/mixedEvent", axisChargePairType, axisInclusiveDeltaEta, axisInclusiveDeltaPhi); + registry.add("inclusive/event/multiplicity", "Inclusive selected-event multiplicity;#it{N}_{tracks}^{PV};events", HistType::kTH1F, {{200, 0.f, 200.f}}); + registry.add("inclusive/normalization/sameEvent", "Inclusive same-event normalization;counter;entries", HistType::kTH1F, {axisCounter}); + registry.add("inclusive/normalization/mixedEvent", "Inclusive mixed-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); + + registry.add("normalization/sameEvent", "Same-event normalization counters;counter;weighted entries", HistType::kTH1F, {axisCounter}); + registry.add("normalization/mixedEventNominal", "Nominal mixed-event normalization counters;counter;weighted entries", HistType::kTH1F, {axisCounter}); + registry.add("normalization/mixedEventLoose", "Loose mixed-event normalization counters;counter;weighted entries", HistType::kTH1F, {axisCounter}); + registry.add("normalization/mixedEventRandom", "Random mixed-event normalization counters;counter;weighted entries", HistType::kTH1F, {axisCounter}); + registry.add("mixedEvent/jetPtDifference", "Mixed-jet matching QA;|#Delta#it{p}_{T,jet}| (GeV/#it{c});jet pairs", HistType::kTH1F, {{80, 0.f, 40.f}}); + registry.add("mixedEvent/jetEtaDifference", "Mixed-jet matching QA;|#Delta#eta_{jet}|;jet pairs", HistType::kTH1F, {{40, 0.f, 1.f}}); + registry.add("mixedEvent/jetNchDifference", "Mixed-jet matching QA;|#Delta#it{N}_{ch}^{jet}|;jet pairs", HistType::kTH1F, {{61, -0.5f, 60.5f}}); + + registry.add("mc/jet/allDetectorPt", "Accepted detector-level charged jets;#it{p}_{T,jet}^{det} (GeV/#it{c});jets", HistType::kTH1F, {axisJetPt}); + registry.add("mc/jet/allParticlePt", "Accepted particle-level charged jets;#it{p}_{T,jet}^{part} (GeV/#it{c});jets", HistType::kTH1F, {axisJetPt}); + registry.add("mc/jet/matchedDetectorPt", "Matched detector-level charged jets;#it{p}_{T,jet}^{det} (GeV/#it{c});jets", HistType::kTH1F, {axisJetPt}); + registry.add("mc/jet/matchedParticlePt", "Matched particle-level charged jets;#it{p}_{T,jet}^{part} (GeV/#it{c});jets", HistType::kTH1F, {axisJetPt}); + registry.add("mc/jet/fakeDetectorPt", "Unmatched detector-level charged jets;#it{p}_{T,jet}^{det} (GeV/#it{c});jets", HistType::kTH1F, {axisJetPt}); + registry.add("mc/jet/responsePt", "Charged-jet response;#it{p}_{T,jet}^{part} (GeV/#it{c});#it{p}_{T,jet}^{det} (GeV/#it{c})", HistType::kTH2F, {axisJetPt, axisJetPt}); + registry.add("mc/jet/relativePtResolution", "Charged-jet momentum resolution;#it{p}_{T,jet}^{part} (GeV/#it{c});(#it{p}_{T}^{det}-#it{p}_{T}^{part})/#it{p}_{T}^{part}", HistType::kTH2F, {axisJetPt, {120, -1.5f, 1.5f}}); + registry.add("mc/jet/matchingDeltaR", "Detector-particle jet matching distance;#Delta#it{R}_{match};matched jets", HistType::kTH1F, {{100, 0.f, 1.f}}); + registry.add("mc/pid/generatedPrimaryProtonPt", "Generated primary protons in accepted particle-level jets;#it{p}_{T}^{part} (GeV/#it{c});particles", HistType::kTH1F, {axisBaryonPt}); + registry.add("mc/pid/reconstructedTrueProtonPt", "Reconstructed primary true protons before PID;#it{p}_{T}^{rec} (GeV/#it{c});tracks", HistType::kTH1F, {axisBaryonPt}); + registry.add("mc/pid/selectedTrueProtonPt", "PID-selected primary true protons;#it{p}_{T}^{rec} (GeV/#it{c});tracks", HistType::kTH1F, {axisBaryonPt}); + registry.add("mc/pid/selectedCandidateOriginVsPt", "Selected proton-candidate origin;#it{p}_{T}^{rec} (GeV/#it{c});origin", HistType::kTH2F, {axisBaryonPt, axisOrigin}); + registry.add("mc/pid/selectedCandidateTruePdgVsPt", "Selected proton-candidate contamination;#it{p}_{T}^{rec} (GeV/#it{c});true PDG code", HistType::kTH2F, {axisBaryonPt, {801, -4005.f, 4005.f}}); + registry.add("mc/pid/ptResponsePrimaryProton", "Primary-proton momentum response;#it{p}_{T}^{part} (GeV/#it{c});#it{p}_{T}^{rec} (GeV/#it{c})", HistType::kTH2F, {axisBaryonPt, axisBaryonPt}); + addPairHistograms("mc/particle/sameEvent", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + addPairHistograms("mc/particle/mixedEvent", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + registry.add("mc/particle/normalizationSameEvent", "Particle-level same-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); + registry.add("mc/particle/normalizationMixedEvent", "Particle-level mixed-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); + + setHistogramLabels(); + } + + void addPairHistograms(std::string const& directory, + AxisSpec const& axisPairType, + AxisSpec const& axisJetPt, + AxisSpec const& axisDeltaEta, + AxisSpec const& axisDeltaPhi, + AxisSpec const& axisDeltaR) + { + registry.add(directory + "/sameSignDeltaEtaDeltaPhi", "Same-sign baryon pairs;#Delta#eta;#Delta#varphi", HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}); + registry.add(directory + "/oppositeSignDeltaEtaDeltaPhi", "Opposite-sign baryon pairs;#Delta#eta;#Delta#varphi", HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}); + registry.add(directory + "/sameSignDeltaR", "Same-sign baryon pairs;#Delta#it{R};weighted pairs", HistType::kTH1F, {axisDeltaR}); + registry.add(directory + "/oppositeSignDeltaR", "Opposite-sign baryon pairs;#Delta#it{R};weighted pairs", HistType::kTH1F, {axisDeltaR}); + registry.add(directory + "/pairSparse", "Baryon-pair distribution;pair type;#it{p}_{T,jet};#Delta#eta;#Delta#varphi;#Delta#it{R}", HistType::kTHnSparseF, {axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR}); + } + + void addInclusivePairHistograms(std::string const& directory, + AxisSpec const& axisChargePairType, + AxisSpec const& axisDeltaEta, + AxisSpec const& axisDeltaPhi) + { + registry.add(directory + "/sameSignDeltaEtaDeltaPhi", "Inclusive same-sign baryon pairs;#Delta#eta;#Delta#varphi", HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}); + registry.add(directory + "/oppositeSignDeltaEtaDeltaPhi", "Inclusive opposite-sign baryon pairs;#Delta#eta;#Delta#varphi", HistType::kTH2F, {axisDeltaEta, axisDeltaPhi}); + registry.add(directory + "/chargePairSparse", "Inclusive charge-separated baryon pairs;charge pair type;#Delta#eta;#Delta#varphi", HistType::kTHnSparseF, {axisChargePairType, axisDeltaEta, axisDeltaPhi}); + } + + void setHistogramLabels() + { + auto eventCutFlow = registry.get(HIST("event/cutFlow")); + eventCutFlow->GetXaxis()->SetBinLabel(1, "all"); + eventCutFlow->GetXaxis()->SetBinLabel(2, "z vertex"); + eventCutFlow->GetXaxis()->SetBinLabel(3, "occupancy"); + eventCutFlow->GetXaxis()->SetBinLabel(4, "event selection"); + + auto setNormalizationLabels = [](std::shared_ptr const& histogram) { + histogram->GetXaxis()->SetBinLabel(1, "events or event pairs"); + histogram->GetXaxis()->SetBinLabel(2, "selected jets or jet pairs"); + histogram->GetXaxis()->SetBinLabel(3, "same-sign pairs"); + histogram->GetXaxis()->SetBinLabel(4, "opposite-sign pairs"); + histogram->GetXaxis()->SetBinLabel(5, "near-side same-sign pairs"); + histogram->GetXaxis()->SetBinLabel(6, "near-side opposite-sign pairs"); + }; + setNormalizationLabels(registry.get(HIST("normalization/sameEvent"))); + setNormalizationLabels(registry.get(HIST("normalization/mixedEventNominal"))); + setNormalizationLabels(registry.get(HIST("normalization/mixedEventLoose"))); + setNormalizationLabels(registry.get(HIST("normalization/mixedEventRandom"))); + } + + template + bool isSelectedCollision(CollisionType const& collision, bool fillCutFlow) + { + if (fillCutFlow) { + registry.fill(HIST("event/cutFlow"), 1.f); + } + if (std::abs(collision.posZ()) > eventCuts.vertexZMax.value) { + return false; + } + if (fillCutFlow) { + registry.fill(HIST("event/cutFlow"), 2.f); + } + if (collision.trackOccupancyInTimeRange() < eventCuts.occupancyMin.value || collision.trackOccupancyInTimeRange() > eventCuts.occupancyMax.value) { + return false; + } + if (fillCutFlow) { + registry.fill(HIST("event/cutFlow"), 3.f); + } + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, true, eventCuts.applyRctSelection.value)) { + return false; + } + if (fillCutFlow) { + registry.fill(HIST("event/cutFlow"), 4.f); + } + return true; + } + + template + bool isSelectedJetGeometry(JetType const& jet) const + { + if (jet.r() != static_cast(std::lround(jetCuts.jetRadius.value * 100.f)) || std::abs(jet.eta()) > jetCuts.jetEtaMax.value) { + return false; + } + if (jetCuts.areaFractionMin.value >= 0.f) { + float const fullArea = o2::constants::math::PI * jetCuts.jetRadius.value * jetCuts.jetRadius.value; + if (jet.area() / fullArea < jetCuts.areaFractionMin.value) { + return false; + } + } + return true; + } + + template + bool isSelectedJet(JetType const& jet) const + { + return isSelectedJetGeometry(jet) && jet.pt() >= jetCuts.jetPtMin.value && jet.pt() < jetCuts.jetPtMax.value; + } + + template + bool isSelectedTrack(TrackType const& track) const + { + if (trackCuts.requirePvContributor.value && !track.isPVContributor()) { + return false; + } + if (!track.hasITS() || !track.hasTPC() || !track.isGlobalTrack()) { + return false; + } + if (track.tpcNClsCrossedRows() < trackCuts.tpcCrossedRowsMin.value || track.tpcChi2NCl() > trackCuts.tpcChi2Max.value || track.itsChi2NCl() > trackCuts.itsChi2Max.value) { + return false; + } + if (std::abs(track.eta()) > trackCuts.protonEtaMax.value || track.pt() < trackCuts.protonPtMin.value || track.pt() >= trackCuts.protonPtMax.value) { + return false; + } + if (std::abs(track.dcaXY()) > trackCuts.dcaXYMax.value || std::abs(track.dcaZ()) > trackCuts.dcaZMax.value) { + return false; + } + return true; + } + + template + bool isSelectedProton(TrackType const& track) const + { + auto const pidMethod = static_cast(pidCuts.method.value); + float const tpcNSigma = std::abs(track.tpcNSigmaPr()); + bool isProton = false; + + if (pidMethod == ProtonPidMethod::TpcOnly || track.pt() < pidCuts.tofPtThreshold.value) { + isProton = tpcNSigma < pidCuts.tpcNSigmaMax.value; + } else if (track.hasTOF()) { + float const tofNSigma = std::abs(track.tofNSigmaPr()); + if (pidMethod == ProtonPidMethod::TpcTofCircular) { + isProton = std::hypot(tpcNSigma, tofNSigma) < pidCuts.combinedNSigmaMax.value; + } else if (pidMethod == ProtonPidMethod::TpcTofRectangular) { + isProton = tpcNSigma < pidCuts.tpcNSigmaMax.value && tofNSigma < pidCuts.tofNSigmaMax.value; + } + } + + if (!isProton || !pidCuts.rejectPionKaon.value) { + return isProton; + } + + bool const isPionCompatible = std::abs(track.tpcNSigmaPi()) < pidCuts.rejectionNSigma.value && + (track.pt() < pidCuts.tofPtThreshold.value || (track.hasTOF() && std::abs(track.tofNSigmaPi()) < pidCuts.rejectionNSigma.value)); + bool const isKaonCompatible = std::abs(track.tpcNSigmaKa()) < pidCuts.rejectionNSigma.value && + (track.pt() < pidCuts.tofPtThreshold.value || (track.hasTOF() && std::abs(track.tofNSigmaKa()) < pidCuts.rejectionNSigma.value)); + return !isPionCompatible && !isKaonCompatible; + } + + template + std::vector collectBaryons(JetType const& jet, bool fillQa, float weight = 1.f) + { + std::vector baryons; + for (auto const& jetTrack : jet.template tracks_as()) { + auto track = jetTrack.template track_as(); + if (!isSelectedTrack(track)) { + continue; + } + + if (fillQa) { + float const signedMomentum = static_cast(track.sign()) * track.p(); + registry.fill(HIST("pid/before/tpcSignalVsSignedP"), signedMomentum, track.tpcSignal(), weight); + registry.fill(HIST("pid/before/tpcNSigmaPiVsPt"), track.pt(), track.tpcNSigmaPi(), weight); + registry.fill(HIST("pid/before/tpcNSigmaKaVsPt"), track.pt(), track.tpcNSigmaKa(), weight); + registry.fill(HIST("pid/before/tpcNSigmaPrVsPt"), track.pt(), track.tpcNSigmaPr(), weight); + registry.fill(HIST("pid/tofMatchDenominator"), track.pt(), weight); + if (track.hasTOF()) { + registry.fill(HIST("pid/tofMatchNumerator"), track.pt(), weight); + registry.fill(HIST("pid/before/tofBetaVsSignedP"), signedMomentum, track.beta(), weight); + registry.fill(HIST("pid/before/tofNSigmaPiVsPt"), track.pt(), track.tofNSigmaPi(), weight); + registry.fill(HIST("pid/before/tofNSigmaKaVsPt"), track.pt(), track.tofNSigmaKa(), weight); + registry.fill(HIST("pid/before/tofNSigmaPrVsPt"), track.pt(), track.tofNSigmaPr(), weight); + } + } + + float const rapidity = RecoDecay::y(std::array{track.px(), track.py(), track.pz()}, o2::constants::physics::MassProton); + if (std::abs(rapidity) > trackCuts.protonRapidityMax.value || !isSelectedProton(track)) { + continue; + } + + baryons.push_back({track.pt(), track.eta(), track.phi(), rapidity, static_cast(track.sign())}); + if (fillQa) { + registry.fill(HIST("proton/pt"), track.pt(), weight); + registry.fill(HIST("proton/eta"), track.eta(), weight); + registry.fill(HIST("proton/rapidity"), rapidity, weight); + registry.fill(HIST("proton/tpcNSigmaVsPt"), track.pt(), track.tpcNSigmaPr(), weight); + if (track.hasTOF()) { + registry.fill(HIST("proton/tofNSigmaVsPt"), track.pt(), track.tofNSigmaPr(), weight); + } + registry.fill(HIST("proton/dcaXYVsPt"), track.pt(), track.dcaXY(), weight); + registry.fill(HIST("proton/dcaZVsPt"), track.pt(), track.dcaZ(), weight); + registry.fill(HIST("pid/selectedPtByCharge"), static_cast(track.sign()), track.pt(), weight); + } + } + return baryons; + } + + static PairType getPairType(SelectedBaryon const& baryon1, SelectedBaryon const& baryon2) + { + return baryon1.sign == baryon2.sign ? PairType::SameSign : PairType::OppositeSign; + } + + static ChargePairType getChargePairType(SelectedBaryon const& baryon1, SelectedBaryon const& baryon2) + { + if (baryon1.sign > 0 && baryon2.sign > 0) { + return ChargePairType::ProtonProton; + } + if (baryon1.sign < 0 && baryon2.sign < 0) { + return ChargePairType::AntiprotonAntiproton; + } + return ChargePairType::ProtonAntiproton; + } + + template + std::vector collectInclusiveBaryons(TracksType const& tracks) + { + std::vector baryons; + for (auto const& track : tracks) { + if (!isSelectedTrack(track)) { + continue; + } + float const rapidity = RecoDecay::y(std::array{track.px(), track.py(), track.pz()}, o2::constants::physics::MassProton); + if (std::abs(rapidity) > trackCuts.protonRapidityMax.value || !isSelectedProton(track)) { + continue; + } + baryons.push_back({track.pt(), track.eta(), track.phi(), rapidity, static_cast(track.sign())}); + } + return baryons; + } + + template + bool isSelectedInclusiveCollision(CollisionType const& collision) const + { + return std::abs(collision.posZ()) <= eventCuts.vertexZMax.value && collision.sel8() && + collision.selection_bit(aod::evsel::kNoITSROFrameBorder) && + collision.selection_bit(aod::evsel::kNoTimeFrameBorder) && + collision.selection_bit(aod::evsel::kNoSameBunchPileup) && + collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV); + } + + void fillInclusivePair(bool isMixed, + SelectedBaryon const& baryon1, + SelectedBaryon const& baryon2, + float weight) + { + PairType const pairType = getPairType(baryon1, baryon2); + ChargePairType const chargePairType = getChargePairType(baryon1, baryon2); + float const deltaEta = baryon1.eta - baryon2.eta; + float const deltaPhi = RecoDecay::constrainAngle(baryon1.phi - baryon2.phi, -o2::constants::math::PIHalf); + float const reflectedDeltaPhi = RecoDecay::constrainAngle(-deltaPhi, -o2::constants::math::PIHalf); + if (isMixed) { + if (pairType == PairType::SameSign) { + registry.fill(HIST("inclusive/mixedEvent/sameSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("inclusive/mixedEvent/sameSignDeltaEtaDeltaPhi"), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + } else { + registry.fill(HIST("inclusive/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("inclusive/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + } + registry.fill(HIST("inclusive/mixedEvent/chargePairSparse"), static_cast(chargePairType), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("inclusive/mixedEvent/chargePairSparse"), static_cast(chargePairType), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + } else { + if (pairType == PairType::SameSign) { + registry.fill(HIST("inclusive/sameEvent/sameSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("inclusive/sameEvent/sameSignDeltaEtaDeltaPhi"), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + } else { + registry.fill(HIST("inclusive/sameEvent/oppositeSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("inclusive/sameEvent/oppositeSignDeltaEtaDeltaPhi"), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + } + registry.fill(HIST("inclusive/sameEvent/chargePairSparse"), static_cast(chargePairType), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("inclusive/sameEvent/chargePairSparse"), static_cast(chargePairType), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + } + } + + template + bool isMatchedJetPair(float jetPtDifference, int jetNchDifference, float jetEtaDifference) const + { + if constexpr (Reference == MixedReference::Random) { + return true; + } else if constexpr (Reference == MixedReference::Nominal) { + return jetPtDifference <= mixing.nominalJetPtDifferenceMax.value && + jetNchDifference <= mixing.nominalJetNchDifferenceMax.value && + jetEtaDifference <= mixing.nominalJetEtaDifferenceMax.value; + } else { + return jetPtDifference <= mixing.looseJetPtDifferenceMax.value && + jetNchDifference <= mixing.looseJetNchDifferenceMax.value && + jetEtaDifference <= mixing.looseJetEtaDifferenceMax.value; + } + } + + template + void fillPair(PairType pairType, ChargePairType chargePairType, float jetPt, float deltaEta, float deltaPhi, float deltaR, float weight) + { + float const pairTypeValue = static_cast(pairType); + auto fillSymmetrized = [&](auto const& sameSignName, auto const& oppositeSignName, auto const& sparseName) { + if (pairType == PairType::SameSign) { + registry.fill(sameSignName, deltaEta, deltaPhi, 0.5f * weight); + registry.fill(sameSignName, -deltaEta, -deltaPhi, 0.5f * weight); + } else { + registry.fill(oppositeSignName, deltaEta, deltaPhi, 0.5f * weight); + registry.fill(oppositeSignName, -deltaEta, -deltaPhi, 0.5f * weight); + } + registry.fill(sparseName, pairTypeValue, jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(sparseName, pairTypeValue, jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + }; + + if constexpr (Reference == MixedReference::Nominal) { + fillSymmetrized(HIST("mixedEvent/nominal/sameSignDeltaEtaDeltaPhi"), HIST("mixedEvent/nominal/oppositeSignDeltaEtaDeltaPhi"), HIST("mixedEvent/nominal/pairSparse")); + if (pairType == PairType::SameSign) { + registry.fill(HIST("mixedEvent/nominal/sameSignDeltaR"), deltaR, weight); + } else { + registry.fill(HIST("mixedEvent/nominal/oppositeSignDeltaR"), deltaR, weight); + } + registry.fill(HIST("mixedEvent/nominal/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("mixedEvent/nominal/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } else if constexpr (Reference == MixedReference::Loose) { + fillSymmetrized(HIST("mixedEvent/loose/sameSignDeltaEtaDeltaPhi"), HIST("mixedEvent/loose/oppositeSignDeltaEtaDeltaPhi"), HIST("mixedEvent/loose/pairSparse")); + if (pairType == PairType::SameSign) { + registry.fill(HIST("mixedEvent/loose/sameSignDeltaR"), deltaR, weight); + } else { + registry.fill(HIST("mixedEvent/loose/oppositeSignDeltaR"), deltaR, weight); + } + registry.fill(HIST("mixedEvent/loose/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("mixedEvent/loose/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } else { + fillSymmetrized(HIST("mixedEvent/random/sameSignDeltaEtaDeltaPhi"), HIST("mixedEvent/random/oppositeSignDeltaEtaDeltaPhi"), HIST("mixedEvent/random/pairSparse")); + if (pairType == PairType::SameSign) { + registry.fill(HIST("mixedEvent/random/sameSignDeltaR"), deltaR, weight); + } else { + registry.fill(HIST("mixedEvent/random/oppositeSignDeltaR"), deltaR, weight); + } + registry.fill(HIST("mixedEvent/random/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("mixedEvent/random/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } + } + + template + void fillMixedNormalization(float counter, float weight) + { + if constexpr (Reference == MixedReference::Nominal) { + registry.fill(HIST("normalization/mixedEventNominal"), counter, weight); + } else if constexpr (Reference == MixedReference::Loose) { + registry.fill(HIST("normalization/mixedEventLoose"), counter, weight); + } else { + registry.fill(HIST("normalization/mixedEventRandom"), counter, weight); + } + } + + void fillSameEventPair(PairType pairType, ChargePairType chargePairType, float jetPt, float deltaEta, float deltaPhi, float deltaR, float weight) + { + float const pairTypeValue = static_cast(pairType); + if (pairType == PairType::SameSign) { + registry.fill(HIST("sameEvent/sameSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("sameEvent/sameSignDeltaEtaDeltaPhi"), -deltaEta, -deltaPhi, 0.5f * weight); + registry.fill(HIST("sameEvent/sameSignDeltaR"), deltaR, weight); + } else { + registry.fill(HIST("sameEvent/oppositeSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("sameEvent/oppositeSignDeltaEtaDeltaPhi"), -deltaEta, -deltaPhi, 0.5f * weight); + registry.fill(HIST("sameEvent/oppositeSignDeltaR"), deltaR, weight); + } + registry.fill(HIST("sameEvent/pairSparse"), pairTypeValue, jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("sameEvent/pairSparse"), pairTypeValue, jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("sameEvent/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("sameEvent/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } + + template + std::vector collectParticleBaryons(JetType const& jet, bool fillQa, float weight) + { + std::vector baryons; + for (auto const& particle : jet.template tracks_as()) { + if (std::abs(particle.pdgCode()) != PDG_t::kProton || !particle.isPhysicalPrimary()) { + continue; + } + if (particle.pt() < trackCuts.protonPtMin.value || particle.pt() >= trackCuts.protonPtMax.value || + std::abs(particle.eta()) > trackCuts.protonEtaMax.value || std::abs(particle.y()) > trackCuts.protonRapidityMax.value) { + continue; + } + int8_t const sign = particle.pdgCode() > 0 ? 1 : -1; + baryons.push_back({particle.pt(), particle.eta(), particle.phi(), particle.y(), sign}); + if (fillQa) { + registry.fill(HIST("mc/pid/generatedPrimaryProtonPt"), particle.pt(), weight); + } + } + return baryons; + } + + void fillParticlePair(bool isMixed, PairType pairType, float jetPt, float deltaEta, float deltaPhi, float deltaR, float weight) + { + if (isMixed) { + if (pairType == PairType::SameSign) { + registry.fill(HIST("mc/particle/mixedEvent/sameSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("mc/particle/mixedEvent/sameSignDeltaEtaDeltaPhi"), -deltaEta, -deltaPhi, 0.5f * weight); + registry.fill(HIST("mc/particle/mixedEvent/sameSignDeltaR"), deltaR, weight); + } else { + registry.fill(HIST("mc/particle/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("mc/particle/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), -deltaEta, -deltaPhi, 0.5f * weight); + registry.fill(HIST("mc/particle/mixedEvent/oppositeSignDeltaR"), deltaR, weight); + } + registry.fill(HIST("mc/particle/mixedEvent/pairSparse"), static_cast(pairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("mc/particle/mixedEvent/pairSparse"), static_cast(pairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } else { + if (pairType == PairType::SameSign) { + registry.fill(HIST("mc/particle/sameEvent/sameSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("mc/particle/sameEvent/sameSignDeltaEtaDeltaPhi"), -deltaEta, -deltaPhi, 0.5f * weight); + registry.fill(HIST("mc/particle/sameEvent/sameSignDeltaR"), deltaR, weight); + } else { + registry.fill(HIST("mc/particle/sameEvent/oppositeSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("mc/particle/sameEvent/oppositeSignDeltaEtaDeltaPhi"), -deltaEta, -deltaPhi, 0.5f * weight); + registry.fill(HIST("mc/particle/sameEvent/oppositeSignDeltaR"), deltaR, weight); + } + registry.fill(HIST("mc/particle/sameEvent/pairSparse"), static_cast(pairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("mc/particle/sameEvent/pairSparse"), static_cast(pairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } + } + + template + void fillMCPidQa(JetType const& jet, float weight) + { + for (auto const& jetTrack : jet.template tracks_as()) { + auto track = jetTrack.template track_as(); + if (!isSelectedTrack(track) || !track.has_mcParticle()) { + continue; + } + auto particle = track.template mcParticle_as(); + bool const isTruePrimaryProton = std::abs(particle.pdgCode()) == PDG_t::kProton && particle.isPhysicalPrimary(); + if (isTruePrimaryProton) { + registry.fill(HIST("mc/pid/reconstructedTrueProtonPt"), track.pt(), weight); + } + float const rapidity = RecoDecay::y(std::array{track.px(), track.py(), track.pz()}, o2::constants::physics::MassProton); + if (std::abs(rapidity) > trackCuts.protonRapidityMax.value || !isSelectedProton(track)) { + continue; + } + int origin = 2; + if (particle.isPhysicalPrimary()) { + origin = 1; + } else if (particle.getProcess() == TMCProcess::kPDecay) { + origin = 3; + } + registry.fill(HIST("mc/pid/selectedCandidateOriginVsPt"), track.pt(), static_cast(origin), weight); + registry.fill(HIST("mc/pid/selectedCandidateTruePdgVsPt"), track.pt(), static_cast(particle.pdgCode()), weight); + if (isTruePrimaryProton) { + registry.fill(HIST("mc/pid/selectedTrueProtonPt"), track.pt(), weight); + registry.fill(HIST("mc/pid/ptResponsePrimaryProton"), particle.pt(), track.pt(), weight); + } + } + } + + template + void fillReconstructedSameEvent(CollisionType const& collision, JetsType const& jets, float weight = 1.f) + { + if (!isSelectedCollision(collision, true)) { + return; + } + + registry.fill(HIST("event/zVertex"), collision.posZ(), weight); + registry.fill(HIST("event/multiplicityFT0C"), collision.multFT0C(), weight); + registry.fill(HIST("normalization/sameEvent"), 1.f, weight); + + for (auto const& jet : jets) { + if (isSelectedJetGeometry(jet)) { + registry.fill(HIST("jet/ptBeforePtSelection"), jet.pt(), weight); + } + if (!isSelectedJet(jet)) { + continue; + } + + registry.fill(HIST("normalization/sameEvent"), 2.f, weight); + registry.fill(HIST("jet/pt"), jet.pt(), weight); + registry.fill(HIST("jet/eta"), jet.eta(), weight); + registry.fill(HIST("jet/phi"), jet.phi(), weight); + registry.fill(HIST("jet/nConstituents"), static_cast(jet.tracksIds().size()), weight); + registry.fill(HIST("jet/nConstituentsVsPt"), jet.pt(), static_cast(jet.tracksIds().size()), weight); + + auto baryons = collectBaryons(jet, true, weight); + int numberOfProtons = 0; + int numberOfAntiprotons = 0; + for (auto const& baryon : baryons) { + baryon.sign > 0 ? ++numberOfProtons : ++numberOfAntiprotons; + } + registry.fill(HIST("jet/nSelectedBaryons"), static_cast(baryons.size()), weight); + registry.fill(HIST("jet/nProtons"), numberOfProtons, weight); + registry.fill(HIST("jet/nAntiprotons"), numberOfAntiprotons, weight); + registry.fill(HIST("jet/nProtonsVsAntiprotons"), numberOfProtons, numberOfAntiprotons, weight); + registry.fill(HIST("jet/nBaryonsVsPt"), jet.pt(), static_cast(baryons.size()), weight); + + for (size_t first = 0; first < baryons.size(); ++first) { + for (size_t second = first + 1; second < baryons.size(); ++second) { + PairType const pairType = getPairType(baryons[first], baryons[second]); + ChargePairType const chargePairType = getChargePairType(baryons[first], baryons[second]); + float const deltaEta = baryons[first].eta - baryons[second].eta; + float const deltaPhi = RecoDecay::constrainAngle(baryons[first].phi - baryons[second].phi, -o2::constants::math::PI); + float const deltaR = std::hypot(deltaEta, deltaPhi); + fillSameEventPair(pairType, chargePairType, jet.pt(), deltaEta, deltaPhi, deltaR, weight); + if (pairType == PairType::SameSign) { + registry.fill(HIST("jet/sameSignPairsVsPt"), jet.pt(), weight); + } else { + registry.fill(HIST("jet/oppositeSignPairsVsPt"), jet.pt(), weight); + } + + float const counterBin = pairType == PairType::SameSign ? 3.f : 4.f; + registry.fill(HIST("normalization/sameEvent"), counterBin, weight); + if (std::abs(deltaPhi) < axes.nearSideDeltaPhiMax.value && std::abs(deltaEta) < axes.nearSideDeltaEtaMax.value) { + registry.fill(HIST("normalization/sameEvent"), pairType == PairType::SameSign ? 5.f : 6.f, weight); + } + } + } + } + } + + void processSameEvent(aod::JetCollision const& collision, + ChargedJetsWithConstituents const& jets, + JetTracksWithOriginal const&, + BaryonTracks const&) + { + fillReconstructedSameEvent(collision, jets); + } + PROCESS_SWITCH(BaryonPairInJets, processSameEvent, "Same-event identified baryon pairs in charged jets", true); + + template + void fillMixedJetPair(JetType1 const& jet1, JetType2 const& jet2, float weight) + { + float const jetPtDifference = std::abs(jet1.pt() - jet2.pt()); + int const jetNchDifference = std::abs(static_cast(jet1.tracksIds().size()) - static_cast(jet2.tracksIds().size())); + float const jetEtaDifference = std::abs(jet1.eta() - jet2.eta()); + if (!isMatchedJetPair(jetPtDifference, jetNchDifference, jetEtaDifference)) { + return; + } + + if constexpr (Reference == MixedReference::Random) { + registry.fill(HIST("mixedEvent/jetPtDifference"), jetPtDifference, weight); + registry.fill(HIST("mixedEvent/jetEtaDifference"), jetEtaDifference, weight); + registry.fill(HIST("mixedEvent/jetNchDifference"), jetNchDifference, weight); + } + + auto baryons1 = collectBaryons(jet1, false); + auto baryons2 = collectBaryons(jet2, false); + float const mixedJetPt = 0.5f * (jet1.pt() + jet2.pt()); + fillMixedNormalization(2.f, weight); + + for (auto const& baryon1 : baryons1) { + float const localEta1 = baryon1.eta - jet1.eta(); + float const localPhi1 = RecoDecay::constrainAngle(baryon1.phi - jet1.phi(), -o2::constants::math::PI); + for (auto const& baryon2 : baryons2) { + float const localEta2 = baryon2.eta - jet2.eta(); + float const localPhi2 = RecoDecay::constrainAngle(baryon2.phi - jet2.phi(), -o2::constants::math::PI); + float const deltaEta = localEta1 - localEta2; + float const deltaPhi = RecoDecay::constrainAngle(localPhi1 - localPhi2, -o2::constants::math::PI); + float const deltaR = std::hypot(deltaEta, deltaPhi); + PairType const pairType = getPairType(baryon1, baryon2); + ChargePairType const chargePairType = getChargePairType(baryon1, baryon2); + fillPair(pairType, chargePairType, mixedJetPt, deltaEta, deltaPhi, deltaR, weight); + fillMixedNormalization(pairType == PairType::SameSign ? 3.f : 4.f, weight); + if (std::abs(deltaPhi) < axes.nearSideDeltaPhiMax.value && std::abs(deltaEta) < axes.nearSideDeltaEtaMax.value) { + fillMixedNormalization(pairType == PairType::SameSign ? 5.f : 6.f, weight); + } + } + } + } + + void processMixedEvent(aod::JetCollisions const& collisions, + ChargedJetsWithConstituents const& jets, + JetTracksWithOriginal const&, + BaryonTracks const&) + { + MixingBinning mixingBinning{{mixing.vertexZBins, mixing.multiplicityBins}}; + auto jetsTuple = std::make_tuple(jets); + SameKindPair mixedCollisionPairs{ + mixingBinning, mixing.numberOfMixedEvents.value, -1, collisions, jetsTuple, &mixingCache}; + + for (auto iterator = mixedCollisionPairs.begin(); iterator != mixedCollisionPairs.end(); ++iterator) { + auto& [collision1, jets1, collision2, jets2] = *iterator; + if (!isSelectedCollision(collision1, false) || !isSelectedCollision(collision2, false)) { + continue; + } + + float const weight = 1.f / static_cast(iterator.currentWindowNeighbours()); + registry.fill(HIST("normalization/mixedEventNominal"), 1.f, weight); + registry.fill(HIST("normalization/mixedEventLoose"), 1.f, weight); + registry.fill(HIST("normalization/mixedEventRandom"), 1.f, weight); + + for (auto const& jet1 : jets1) { + if (!isSelectedJet(jet1)) { + continue; + } + for (auto const& jet2 : jets2) { + if (!isSelectedJet(jet2)) { + continue; + } + fillMixedJetPair(jet1, jet2, weight); + fillMixedJetPair(jet1, jet2, weight); + fillMixedJetPair(jet1, jet2, weight); + } + } + } + } + PROCESS_SWITCH(BaryonPairInJets, processMixedEvent, "Mixed-event identified baryon pairs in charged jets", true); + + void processMCD(aod::JetCollisionsMCD::iterator const& collision, + aod::JetMcCollisions const&, + ChargedMCDJetsWithConstituents const& jets, + JetTracksMCDWithOriginal const&, + BaryonTracksMC const&, + aod::McParticles const&) + { + float const weight = collision.mcCollision().weight(); + fillReconstructedSameEvent(collision, jets, weight); + if (!isSelectedCollision(collision, false)) { + return; + } + for (auto const& jet : jets) { + if (isSelectedJet(jet)) { + fillMCPidQa(jet, weight); + } + } + } + PROCESS_SWITCH(BaryonPairInJets, processMCD, "Detector-level MC reconstructed baryon pairs and PID truth QA", false); + + void processMixedEventMCD(aod::JetCollisionsMCD const& collisions, + aod::JetMcCollisions const&, + ChargedMCDJetsWithConstituents const& jets, + JetTracksMCDWithOriginal const&, + BaryonTracksMC const&, + aod::McParticles const&) + { + MixingBinning mixingBinning{{mixing.vertexZBins, mixing.multiplicityBins}}; + auto jetsTuple = std::make_tuple(jets); + SameKindPair mixedCollisionPairs{ + mixingBinning, mixing.numberOfMixedEvents.value, -1, collisions, jetsTuple, &mcdMixingCache}; + + for (auto iterator = mixedCollisionPairs.begin(); iterator != mixedCollisionPairs.end(); ++iterator) { + auto& [collision1, jets1, collision2, jets2] = *iterator; + if (!isSelectedCollision(collision1, false) || !isSelectedCollision(collision2, false)) { + continue; + } + float const weight = collision1.mcCollision().weight() / static_cast(iterator.currentWindowNeighbours()); + registry.fill(HIST("normalization/mixedEventNominal"), 1.f, weight); + registry.fill(HIST("normalization/mixedEventLoose"), 1.f, weight); + registry.fill(HIST("normalization/mixedEventRandom"), 1.f, weight); + for (auto const& jet1 : jets1) { + if (!isSelectedJet(jet1)) { + continue; + } + for (auto const& jet2 : jets2) { + if (!isSelectedJet(jet2)) { + continue; + } + fillMixedJetPair(jet1, jet2, weight); + fillMixedJetPair(jet1, jet2, weight); + fillMixedJetPair(jet1, jet2, weight); + } + } + } + } + PROCESS_SWITCH(BaryonPairInJets, processMixedEventMCD, "Detector-level MC mixed-event baryon pairs", false); + + void processMCP(aod::JetMcCollision const& collision, + ChargedMCPJetsWithConstituents const& jets, + aod::JetParticles const&) + { + if (std::abs(collision.posZ()) > eventCuts.vertexZMax.value) { + return; + } + float const weight = collision.weight(); + registry.fill(HIST("mc/particle/normalizationSameEvent"), 1.f, weight); + for (auto const& jet : jets) { + if (!isSelectedJet(jet)) { + continue; + } + registry.fill(HIST("mc/jet/allParticlePt"), jet.pt(), weight); + registry.fill(HIST("mc/particle/normalizationSameEvent"), 2.f, weight); + auto baryons = collectParticleBaryons(jet, true, weight); + for (size_t first = 0; first < baryons.size(); ++first) { + for (size_t second = first + 1; second < baryons.size(); ++second) { + PairType const pairType = getPairType(baryons[first], baryons[second]); + float const deltaEta = baryons[first].eta - baryons[second].eta; + float const deltaPhi = RecoDecay::constrainAngle(baryons[first].phi - baryons[second].phi, -o2::constants::math::PI); + float const deltaR = std::hypot(deltaEta, deltaPhi); + fillParticlePair(false, pairType, jet.pt(), deltaEta, deltaPhi, deltaR, weight); + registry.fill(HIST("mc/particle/normalizationSameEvent"), pairType == PairType::SameSign ? 3.f : 4.f, weight); + } + } + } + } + PROCESS_SWITCH(BaryonPairInJets, processMCP, "Particle-level MC baryon pairs", false); + + void processMixedEventMCP(aod::JetMcCollisions const& collisions, + ChargedMCPJetsWithConstituents const& jets, + aod::JetParticles const&) + { + MCPMixingBinning mixingBinning{{mixing.vertexZBins, mixing.multiplicityBins}}; + auto jetsTuple = std::make_tuple(jets); + SameKindPair mixedCollisionPairs{ + mixingBinning, mixing.numberOfMixedEvents.value, -1, collisions, jetsTuple, &mcpMixingCache}; + for (auto iterator = mixedCollisionPairs.begin(); iterator != mixedCollisionPairs.end(); ++iterator) { + auto& [collision1, jets1, collision2, jets2] = *iterator; + if (std::abs(collision1.posZ()) > eventCuts.vertexZMax.value || std::abs(collision2.posZ()) > eventCuts.vertexZMax.value) { + continue; + } + float const weight = collision1.weight() / static_cast(iterator.currentWindowNeighbours()); + registry.fill(HIST("mc/particle/normalizationMixedEvent"), 1.f, weight); + for (auto const& jet1 : jets1) { + if (!isSelectedJet(jet1)) { + continue; + } + for (auto const& jet2 : jets2) { + if (!isSelectedJet(jet2)) { + continue; + } + auto baryons1 = collectParticleBaryons(jet1, false, weight); + auto baryons2 = collectParticleBaryons(jet2, false, weight); + registry.fill(HIST("mc/particle/normalizationMixedEvent"), 2.f, weight); + for (auto const& baryon1 : baryons1) { + float const localEta1 = baryon1.eta - jet1.eta(); + float const localPhi1 = RecoDecay::constrainAngle(baryon1.phi - jet1.phi(), -o2::constants::math::PI); + for (auto const& baryon2 : baryons2) { + float const localEta2 = baryon2.eta - jet2.eta(); + float const localPhi2 = RecoDecay::constrainAngle(baryon2.phi - jet2.phi(), -o2::constants::math::PI); + float const deltaEta = localEta1 - localEta2; + float const deltaPhi = RecoDecay::constrainAngle(localPhi1 - localPhi2, -o2::constants::math::PI); + float const deltaR = std::hypot(deltaEta, deltaPhi); + PairType const pairType = getPairType(baryon1, baryon2); + fillParticlePair(true, pairType, 0.5f * (jet1.pt() + jet2.pt()), deltaEta, deltaPhi, deltaR, weight); + registry.fill(HIST("mc/particle/normalizationMixedEvent"), pairType == PairType::SameSign ? 3.f : 4.f, weight); + } + } + } + } + } + } + PROCESS_SWITCH(BaryonPairInJets, processMixedEventMCP, "Particle-level MC mixed-event baryon pairs", false); + + void processMatched(aod::JetCollisionsMCD::iterator const& collision, + aod::JetMcCollisions const&, + MatchedChargedMCDJets const& detectorJets, + MatchedChargedMCPJets const&, + JetTracksMCDWithOriginal const&, + aod::JetParticles const&) + { + if (!isSelectedCollision(collision, false)) { + return; + } + float const weight = collision.mcCollision().weight(); + for (auto const& detectorJet : detectorJets) { + if (!isSelectedJet(detectorJet)) { + continue; + } + registry.fill(HIST("mc/jet/allDetectorPt"), detectorJet.pt(), weight); + bool hasAcceptedMatch = false; + for (auto const& particleJet : detectorJet.template matchedJetGeo_as()) { + if (!isSelectedJetGeometry(particleJet)) { + continue; + } + hasAcceptedMatch = true; + registry.fill(HIST("mc/jet/matchedDetectorPt"), detectorJet.pt(), weight); + registry.fill(HIST("mc/jet/matchedParticlePt"), particleJet.pt(), weight); + registry.fill(HIST("mc/jet/responsePt"), particleJet.pt(), detectorJet.pt(), weight); + if (particleJet.pt() > 0.f) { + registry.fill(HIST("mc/jet/relativePtResolution"), particleJet.pt(), (detectorJet.pt() - particleJet.pt()) / particleJet.pt(), weight); + } + float const deltaEta = detectorJet.eta() - particleJet.eta(); + float const deltaPhi = RecoDecay::constrainAngle(detectorJet.phi() - particleJet.phi(), -o2::constants::math::PI); + registry.fill(HIST("mc/jet/matchingDeltaR"), std::hypot(deltaEta, deltaPhi), weight); + break; + } + if (!hasAcceptedMatch) { + registry.fill(HIST("mc/jet/fakeDetectorPt"), detectorJet.pt(), weight); + } + } + } + PROCESS_SWITCH(BaryonPairInJets, processMatched, "Matched detector-particle charged-jet response", false); + + void processInclusiveSameEvent(InclusiveCollisions::iterator const& collision, BaryonTracks const& tracks) + { + if (!isSelectedInclusiveCollision(collision)) { + return; + } + registry.fill(HIST("inclusive/event/multiplicity"), collision.multNTracksPV()); + registry.fill(HIST("inclusive/normalization/sameEvent"), 1.f); + auto baryons = collectInclusiveBaryons(tracks); + registry.fill(HIST("inclusive/normalization/sameEvent"), 2.f, static_cast(baryons.size())); + for (size_t first = 0; first < baryons.size(); ++first) { + for (size_t second = first + 1; second < baryons.size(); ++second) { + fillInclusivePair(false, baryons[first], baryons[second], 1.f); + registry.fill(HIST("inclusive/normalization/sameEvent"), getPairType(baryons[first], baryons[second]) == PairType::SameSign ? 3.f : 4.f); + } + } + } + PROCESS_SWITCH(BaryonPairInJets, processInclusiveSameEvent, "Inclusive same-event identified baryon pairs", true); + + void processInclusiveMixedEvent(InclusiveCollisions const& collisions, BaryonTracks const& tracks) + { + InclusiveMixingBinning mixingBinning{{mixing.vertexZBins, mixing.multiplicityBins}}; + auto tracksTuple = std::make_tuple(tracks); + SameKindPair mixedCollisionPairs{ + mixingBinning, mixing.numberOfMixedEvents.value, -1, collisions, tracksTuple, &inclusiveMixingCache}; + + for (auto iterator = mixedCollisionPairs.begin(); iterator != mixedCollisionPairs.end(); ++iterator) { + auto& [collision1, tracks1, collision2, tracks2] = *iterator; + if (!isSelectedInclusiveCollision(collision1) || !isSelectedInclusiveCollision(collision2)) { + continue; + } + float const weight = 1.f / static_cast(iterator.currentWindowNeighbours()); + registry.fill(HIST("inclusive/normalization/mixedEvent"), 1.f, weight); + auto baryons1 = collectInclusiveBaryons(tracks1); + auto baryons2 = collectInclusiveBaryons(tracks2); + registry.fill(HIST("inclusive/normalization/mixedEvent"), 2.f, weight * static_cast(baryons1.size() + baryons2.size())); + for (auto const& baryon1 : baryons1) { + for (auto const& baryon2 : baryons2) { + fillInclusivePair(true, baryon1, baryon2, weight); + registry.fill(HIST("inclusive/normalization/mixedEvent"), getPairType(baryon1, baryon2) == PairType::SameSign ? 3.f : 4.f, weight); + } + } + } + } + PROCESS_SWITCH(BaryonPairInJets, processInclusiveMixedEvent, "Inclusive mixed-event identified baryon pairs", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& configContext) +{ + return WorkflowSpec{adaptAnalysisTask(configContext)}; +} From ae88ad993aadeffc18095f3d65cc15472a07c3e4 Mon Sep 17 00:00:00 2001 From: HANSEO PARK Date: Mon, 12 Oct 2026 20:21:38 +0900 Subject: [PATCH 2/2] fix megalinter --- PWGJE/Tasks/baryonPairInJets.cxx | 728 +++++++++++++++++++++++++------ 1 file changed, 586 insertions(+), 142 deletions(-) diff --git a/PWGJE/Tasks/baryonPairInJets.cxx b/PWGJE/Tasks/baryonPairInJets.cxx index bddcf690ae7..80539c2bbfb 100644 --- a/PWGJE/Tasks/baryonPairInJets.cxx +++ b/PWGJE/Tasks/baryonPairInJets.cxx @@ -17,6 +17,7 @@ #include "PWGJE/DataModel/Jet.h" #include "PWGJE/DataModel/JetReducedData.h" +#include "Common/CCDB/RCTSelectionFlags.h" #include "Common/Core/RecoDecay.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" @@ -79,9 +80,9 @@ struct BaryonPairInJets { }; enum class ChargePairType : uint8_t { - ProtonProton = 0, - AntiprotonAntiproton = 1, - ProtonAntiproton = 2 + PositivePositive = 0, + NegativeNegative = 1, + OppositeSign = 2 }; enum class MixedReference : uint8_t { @@ -96,18 +97,31 @@ struct BaryonPairInJets { TpcTofRectangular = 2 }; - struct SelectedBaryon { + enum class ControlSpecies : uint8_t { + Pion = 0, + Kaon = 1 + }; + + // Species-independent selected-particle representation. Keep the full + // kinematics needed by pair and future jet-fragmentation observables. + struct SelectedParticle { float pt; float eta; float phi; - float rapidity; int8_t sign; }; + struct SelectedParticles { + std::vector baryons; + std::vector pions; + std::vector kaons; + }; + struct : ConfigurableGroup { Configurable vertexZMax{"vertexZMax", 10.f, "Maximum absolute collision z vertex (cm)"}; Configurable eventSelections{"eventSelections", "sel8+NoITSROFrameBorder+NoTimeFrameBorder+NoSameBunchPileup+IsGoodZvtxFT0vsPV", "Event selections"}; Configurable applyRctSelection{"applyRctSelection", true, "Apply the PWGJE hadron-PID RCT selection"}; + Configurable rctLabel{"rctLabel", "CBT_hadronPID", "RCT selection label"}; Configurable occupancyMin{"occupancyMin", -999999, "Minimum track occupancy in time range"}; Configurable occupancyMax{"occupancyMax", 999999, "Maximum track occupancy in time range"}; } eventCuts; @@ -123,8 +137,7 @@ struct BaryonPairInJets { struct : ConfigurableGroup { Configurable protonPtMin{"protonPtMin", 0.5f, "Minimum proton transverse momentum (GeV/c)"}; Configurable protonPtMax{"protonPtMax", 2.5f, "Maximum proton transverse momentum (GeV/c)"}; - Configurable protonRapidityMax{"protonRapidityMax", 0.5f, "Maximum absolute proton rapidity"}; - Configurable protonEtaMax{"protonEtaMax", 0.8f, "Maximum absolute proton pseudorapidity before the rapidity cut"}; + Configurable protonEtaMax{"protonEtaMax", 0.8f, "Maximum absolute proton pseudorapidity"}; Configurable requirePvContributor{"requirePvContributor", false, "Require the track to be a PV contributor"}; Configurable tpcCrossedRowsMin{"tpcCrossedRowsMin", 70, "Minimum number of crossed TPC rows"}; Configurable tpcChi2Max{"tpcChi2Max", 4.f, "Maximum TPC chi2 per cluster"}; @@ -141,6 +154,12 @@ struct BaryonPairInJets { Configurable combinedNSigmaMax{"combinedNSigmaMax", 3.f, "Maximum circular TPC-TOF proton n-sigma"}; Configurable rejectPionKaon{"rejectPionKaon", false, "Reject candidates compatible with the pion or kaon hypothesis"}; Configurable rejectionNSigma{"rejectionNSigma", 3.f, "Pion and kaon rejection threshold"}; + Configurable pionTpcNSigmaMax{"pionTpcNSigmaMax", 3.f, "Maximum absolute pion TPC n-sigma"}; + Configurable pionTofNSigmaMax{"pionTofNSigmaMax", 3.f, "Maximum absolute pion TOF n-sigma"}; + Configurable pionCombinedNSigmaMax{"pionCombinedNSigmaMax", 3.f, "Maximum circular TPC-TOF pion n-sigma"}; + Configurable kaonTpcNSigmaMax{"kaonTpcNSigmaMax", 3.f, "Maximum absolute kaon TPC n-sigma"}; + Configurable kaonTofNSigmaMax{"kaonTofNSigmaMax", 3.f, "Maximum absolute kaon TOF n-sigma"}; + Configurable kaonCombinedNSigmaMax{"kaonCombinedNSigmaMax", 3.f, "Maximum circular TPC-TOF kaon n-sigma"}; } pidCuts; struct : ConfigurableGroup { @@ -155,6 +174,14 @@ struct BaryonPairInJets { Configurable looseJetEtaDifferenceMax{"looseJetEtaDifferenceMax", 0.3f, "Maximum jet-eta difference for loose mixing"}; } mixing; + struct : ConfigurableGroup { + Configurable fillQa{"fillQa", true, "Fill event, jet, track, and PID QA histograms"}; + Configurable fillReferenceVariations{"fillReferenceVariations", true, "Fill loose and unmatched mixed-event references"}; + Configurable fillChargeSeparated{"fillChargeSeparated", true, "Fill charge-resolved pair histograms"}; + Configurable fillPionControl{"fillPionControl", false, "Fill same-species pion control correlations"}; + Configurable fillKaonControl{"fillKaonControl", false, "Fill same-species kaon control correlations"}; + } output; + struct : ConfigurableGroup { Configurable nearSideDeltaPhiMax{"nearSideDeltaPhiMax", 0.2f, "Maximum absolute delta phi for near-side pair yields"}; Configurable nearSideDeltaEtaMax{"nearSideDeltaEtaMax", 0.6f, "Maximum absolute delta eta for near-side pair yields"}; @@ -177,18 +204,20 @@ struct BaryonPairInJets { HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; std::vector eventSelectionBits; + o2::aod::rctsel::RCTFlagsChecker inclusiveRctChecker; SliceCache mixingCache; SliceCache mcdMixingCache; SliceCache mcpMixingCache; SliceCache inclusiveMixingCache; using MixingBinning = ColumnBinningPolicy; using MCPMixingBinning = ColumnBinningPolicy; - using InclusiveCollisions = soa::Join; - using InclusiveMixingBinning = ColumnBinningPolicy; + using InclusiveCollisions = soa::Join; + using InclusiveMixingBinning = ColumnBinningPolicy; void init(InitContext const&) { eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventCuts.eventSelections)); + inclusiveRctChecker.init(static_cast(eventCuts.rctLabel)); AxisSpec axisPairType{2, -0.5, 1.5, "pair type"}; AxisSpec axisChargePairType{3, -0.5, 2.5, "charge pair type"}; @@ -211,41 +240,43 @@ struct BaryonPairInJets { AxisSpec axisOrigin{3, 0.5, 3.5, "origin"}; registry.add("event/cutFlow", "Event selection;selection step;entries", HistType::kTH1F, {axisCounter}); - registry.add("event/zVertex", "Selected collision z vertex;#it{z}_{vtx} (cm);entries", HistType::kTH1F, {{80, -20.f, 20.f}}); - registry.add("event/multiplicityFT0C", "Selected collision FT0C multiplicity;FT0C multiplicity;entries", HistType::kTH1F, {{200, 0.f, 2000.f}}); - - registry.add("jet/ptBeforePtSelection", "Charged jets before the jet-pT selection;#it{p}_{T,jet} (GeV/#it{c});entries", HistType::kTH1F, {axisJetPt}); - registry.add("jet/pt", "Selected charged jets;#it{p}_{T,jet} (GeV/#it{c});entries", HistType::kTH1F, {axisJetPt}); - registry.add("jet/eta", "Selected charged jets;#eta_{jet};entries", HistType::kTH1F, {axisJetEta}); - registry.add("jet/phi", "Selected charged jets;#varphi_{jet};entries", HistType::kTH1F, {axisJetPhi}); - registry.add("jet/nConstituents", "Charged-jet constituent multiplicity;#it{N}_{ch}^{jet};entries", HistType::kTH1F, {axisJetNch}); - registry.add("jet/nSelectedBaryons", "Selected baryon multiplicity per jet;#it{N}_{p+#bar{p}};jets", HistType::kTH1F, {axisMultiplicity}); - registry.add("jet/nProtons", "Selected proton multiplicity per jet;#it{N}_{p};jets", HistType::kTH1F, {axisMultiplicity}); - registry.add("jet/nAntiprotons", "Selected antiproton multiplicity per jet;#it{N}_{#bar{p}};jets", HistType::kTH1F, {axisMultiplicity}); - registry.add("jet/nProtonsVsAntiprotons", "Proton and antiproton multiplicities per jet;#it{N}_{p};#it{N}_{#bar{p}}", HistType::kTH2F, {axisMultiplicity, axisMultiplicity}); - registry.add("jet/nConstituentsVsPt", "Constituent multiplicity versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});#it{N}_{ch}^{jet}", HistType::kTH2F, {axisJetPt, axisJetNch}); - registry.add("jet/nBaryonsVsPt", "Selected baryon multiplicity versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});#it{N}_{p+#bar{p}}", HistType::kTH2F, {axisJetPt, axisMultiplicity}); - registry.add("jet/sameSignPairsVsPt", "Same-sign pair count versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});pairs", HistType::kTH1F, {axisJetPt}); - registry.add("jet/oppositeSignPairsVsPt", "Opposite-sign pair count versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});pairs", HistType::kTH1F, {axisJetPt}); - - registry.add("proton/pt", "Selected proton and antiproton candidates;#it{p}_{T} (GeV/#it{c});entries", HistType::kTH1F, {axisBaryonPt}); - registry.add("proton/eta", "Selected proton and antiproton candidates;#eta;entries", HistType::kTH1F, {axisBaryonEta}); - registry.add("proton/rapidity", "Selected proton and antiproton candidates;#it{y};entries", HistType::kTH1F, {axisBaryonRapidity}); - registry.add("proton/tpcNSigmaVsPt", "Proton TPC PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); - registry.add("proton/tofNSigmaVsPt", "Proton TOF PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); - registry.add("proton/dcaXYVsPt", "Selected proton transverse DCA;#it{p}_{T} (GeV/#it{c});DCA_{xy} (cm)", HistType::kTH2F, {axisBaryonPt, axisDca}); - registry.add("proton/dcaZVsPt", "Selected proton longitudinal DCA;#it{p}_{T} (GeV/#it{c});DCA_{z} (cm)", HistType::kTH2F, {axisBaryonPt, axisDca}); - registry.add("pid/before/tpcSignalVsSignedP", "TPC response before proton PID;charge #times #it{p} (GeV/#it{c});TPC signal", HistType::kTH2F, {{120, -6.f, 6.f}, {200, 0.f, 200.f}}); - registry.add("pid/before/tofBetaVsSignedP", "TOF response before proton PID;charge #times #it{p} (GeV/#it{c});#beta", HistType::kTH2F, {{120, -6.f, 6.f}, {120, 0.f, 1.2f}}); - registry.add("pid/before/tpcNSigmaPiVsPt", "Pion TPC PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{#pi}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); - registry.add("pid/before/tpcNSigmaKaVsPt", "Kaon TPC PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{K}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); - registry.add("pid/before/tpcNSigmaPrVsPt", "Proton TPC PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); - registry.add("pid/before/tofNSigmaPiVsPt", "Pion TOF PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{#pi}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); - registry.add("pid/before/tofNSigmaKaVsPt", "Kaon TOF PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{K}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); - registry.add("pid/before/tofNSigmaPrVsPt", "Proton TOF PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); - registry.add("pid/tofMatchDenominator", "Tracks entering PID QA;#it{p}_{T} (GeV/#it{c});tracks", HistType::kTH1F, {axisBaryonPt}); - registry.add("pid/tofMatchNumerator", "Tracks with TOF;#it{p}_{T} (GeV/#it{c});tracks", HistType::kTH1F, {axisBaryonPt}); - registry.add("pid/selectedPtByCharge", "Selected proton candidates;charge;#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{2, -1.5f, 1.5f}, axisBaryonPt}); + if (output.fillQa.value) { + registry.add("event/zVertex", "Selected collision z vertex;#it{z}_{vtx} (cm);entries", HistType::kTH1F, {{80, -20.f, 20.f}}); + registry.add("event/multiplicityFT0C", "Selected collision FT0C multiplicity;FT0C multiplicity;entries", HistType::kTH1F, {{200, 0.f, 2000.f}}); + + registry.add("jet/ptBeforePtSelection", "Charged jets before the jet-pT selection;#it{p}_{T,jet} (GeV/#it{c});entries", HistType::kTH1F, {axisJetPt}); + registry.add("jet/pt", "Selected charged jets;#it{p}_{T,jet} (GeV/#it{c});entries", HistType::kTH1F, {axisJetPt}); + registry.add("jet/eta", "Selected charged jets;#eta_{jet};entries", HistType::kTH1F, {axisJetEta}); + registry.add("jet/phi", "Selected charged jets;#varphi_{jet};entries", HistType::kTH1F, {axisJetPhi}); + registry.add("jet/nConstituents", "Charged-jet constituent multiplicity;#it{N}_{ch}^{jet};entries", HistType::kTH1F, {axisJetNch}); + registry.add("jet/nSelectedBaryons", "Selected baryon multiplicity per jet;#it{N}_{p+#bar{p}};jets", HistType::kTH1F, {axisMultiplicity}); + registry.add("jet/nProtons", "Selected proton multiplicity per jet;#it{N}_{p};jets", HistType::kTH1F, {axisMultiplicity}); + registry.add("jet/nAntiprotons", "Selected antiproton multiplicity per jet;#it{N}_{#bar{p}};jets", HistType::kTH1F, {axisMultiplicity}); + registry.add("jet/nProtonsVsAntiprotons", "Proton and antiproton multiplicities per jet;#it{N}_{p};#it{N}_{#bar{p}}", HistType::kTH2F, {axisMultiplicity, axisMultiplicity}); + registry.add("jet/nConstituentsVsPt", "Constituent multiplicity versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});#it{N}_{ch}^{jet}", HistType::kTH2F, {axisJetPt, axisJetNch}); + registry.add("jet/nBaryonsVsPt", "Selected baryon multiplicity versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});#it{N}_{p+#bar{p}}", HistType::kTH2F, {axisJetPt, axisMultiplicity}); + registry.add("jet/sameSignPairsVsPt", "Same-sign pair count versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});pairs", HistType::kTH1F, {axisJetPt}); + registry.add("jet/oppositeSignPairsVsPt", "Opposite-sign pair count versus jet momentum;#it{p}_{T,jet} (GeV/#it{c});pairs", HistType::kTH1F, {axisJetPt}); + + registry.add("proton/pt", "Selected proton and antiproton candidates;#it{p}_{T} (GeV/#it{c});entries", HistType::kTH1F, {axisBaryonPt}); + registry.add("proton/eta", "Selected proton and antiproton candidates;#eta;entries", HistType::kTH1F, {axisBaryonEta}); + registry.add("proton/rapidity", "Selected proton and antiproton candidates;#it{y};entries", HistType::kTH1F, {axisBaryonRapidity}); + registry.add("proton/tpcNSigmaVsPt", "Proton TPC PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("proton/tofNSigmaVsPt", "Proton TOF PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("proton/dcaXYVsPt", "Selected proton transverse DCA;#it{p}_{T} (GeV/#it{c});DCA_{xy} (cm)", HistType::kTH2F, {axisBaryonPt, axisDca}); + registry.add("proton/dcaZVsPt", "Selected proton longitudinal DCA;#it{p}_{T} (GeV/#it{c});DCA_{z} (cm)", HistType::kTH2F, {axisBaryonPt, axisDca}); + registry.add("pid/before/tpcSignalVsSignedP", "TPC response before proton PID;charge #times #it{p} (GeV/#it{c});TPC signal", HistType::kTH2F, {{120, -6.f, 6.f}, {200, 0.f, 200.f}}); + registry.add("pid/before/tofBetaVsSignedP", "TOF response before proton PID;charge #times #it{p} (GeV/#it{c});#beta", HistType::kTH2F, {{120, -6.f, 6.f}, {120, 0.f, 1.2f}}); + registry.add("pid/before/tpcNSigmaPiVsPt", "Pion TPC PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{#pi}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tpcNSigmaKaVsPt", "Kaon TPC PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{K}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tpcNSigmaPrVsPt", "Proton TPC PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tofNSigmaPiVsPt", "Pion TOF PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{#pi}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tofNSigmaKaVsPt", "Kaon TOF PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{K}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/before/tofNSigmaPrVsPt", "Proton TOF PID before proton PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}^{p}", HistType::kTH2F, {axisBaryonPt, axisNSigma}); + registry.add("pid/tofMatchDenominator", "Tracks entering PID QA;#it{p}_{T} (GeV/#it{c});tracks", HistType::kTH1F, {axisBaryonPt}); + registry.add("pid/tofMatchNumerator", "Tracks with TOF;#it{p}_{T} (GeV/#it{c});tracks", HistType::kTH1F, {axisBaryonPt}); + registry.add("pid/selectedPtByCharge", "Selected proton candidates;charge;#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{2, -1.5f, 1.5f}, axisBaryonPt}); + } addPairHistograms("sameEvent", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); addPairHistograms("mixedEvent/nominal", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); @@ -258,7 +289,7 @@ struct BaryonPairInJets { addInclusivePairHistograms("inclusive/sameEvent", axisChargePairType, axisInclusiveDeltaEta, axisInclusiveDeltaPhi); addInclusivePairHistograms("inclusive/mixedEvent", axisChargePairType, axisInclusiveDeltaEta, axisInclusiveDeltaPhi); - registry.add("inclusive/event/multiplicity", "Inclusive selected-event multiplicity;#it{N}_{tracks}^{PV};events", HistType::kTH1F, {{200, 0.f, 200.f}}); + registry.add("inclusive/event/multiplicity", "Inclusive selected-event multiplicity;FT0C multiplicity;events", HistType::kTH1F, {{200, 0.f, 2000.f}}); registry.add("inclusive/normalization/sameEvent", "Inclusive same-event normalization;counter;entries", HistType::kTH1F, {axisCounter}); registry.add("inclusive/normalization/mixedEvent", "Inclusive mixed-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); @@ -289,9 +320,53 @@ struct BaryonPairInJets { registry.add("mc/particle/normalizationSameEvent", "Particle-level same-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); registry.add("mc/particle/normalizationMixedEvent", "Particle-level mixed-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); + if (output.fillPionControl.value) { + addControlSpeciesHistograms("control/pion", "Pion", axisPairType, axisCounter, axisJetPt, axisBaryonPt, axisNSigma, axisOrigin, axisDeltaEta, axisDeltaPhi, axisDeltaR); + } + if (output.fillKaonControl.value) { + addControlSpeciesHistograms("control/kaon", "Kaon", axisPairType, axisCounter, axisJetPt, axisBaryonPt, axisNSigma, axisOrigin, axisDeltaEta, axisDeltaPhi, axisDeltaR); + } + setHistogramLabels(); } + void addControlSpeciesHistograms(std::string const& directory, + std::string const& speciesName, + AxisSpec const& axisPairType, + AxisSpec const& axisCounter, + AxisSpec const& axisJetPt, + AxisSpec const& axisParticlePt, + AxisSpec const& axisNSigma, + AxisSpec const& axisOrigin, + AxisSpec const& axisDeltaEta, + AxisSpec const& axisDeltaPhi, + AxisSpec const& axisDeltaR) + { + auto addHistogram = [&](std::string const& suffix, std::string const& title, HistType type, std::vector const& histogramAxes) { + registry.add(directory + suffix, title.c_str(), type, histogramAxes); + }; + AxisSpec const axisPdg{801, -4005.f, 4005.f, "true PDG code"}; + addPairHistograms(directory + "/sameEvent", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + addPairHistograms(directory + "/mixedEvent", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + addHistogram("/normalizationSameEvent", speciesName + " same-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); + addHistogram("/normalizationMixedEvent", speciesName + " mixed-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); + if (output.fillQa.value) { + addHistogram("/qa/pt", "Selected " + speciesName + " candidates;#it{p}_{T} (GeV/#it{c});entries", HistType::kTH1F, {axisParticlePt}); + addHistogram("/qa/tpcNSigmaVsPt", speciesName + " TPC PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}", HistType::kTH2F, {axisParticlePt, axisNSigma}); + addHistogram("/qa/tofNSigmaVsPt", speciesName + " TOF PID;#it{p}_{T} (GeV/#it{c});n#sigma_{TOF}", HistType::kTH2F, {axisParticlePt, axisNSigma}); + } + addHistogram("/mc/generatedPrimaryPt", "Generated primary " + speciesName + " candidates;#it{p}_{T}^{part} (GeV/#it{c});particles", HistType::kTH1F, {axisParticlePt}); + addHistogram("/mc/reconstructedTruePt", "Reconstructed true primary " + speciesName + " candidates;#it{p}_{T}^{rec} (GeV/#it{c});tracks", HistType::kTH1F, {axisParticlePt}); + addHistogram("/mc/selectedTruePt", "PID-selected true primary " + speciesName + " candidates;#it{p}_{T}^{rec} (GeV/#it{c});tracks", HistType::kTH1F, {axisParticlePt}); + addHistogram("/mc/selectedCandidateOriginVsPt", "Selected " + speciesName + " candidate origin;#it{p}_{T}^{rec} (GeV/#it{c});origin", HistType::kTH2F, {axisParticlePt, axisOrigin}); + addHistogram("/mc/selectedCandidateTruePdgVsPt", "Selected " + speciesName + " candidate contamination;#it{p}_{T}^{rec} (GeV/#it{c});true PDG code", HistType::kTH2F, {axisParticlePt, axisPdg}); + addHistogram("/mc/ptResponsePrimary", "Primary " + speciesName + " momentum response;#it{p}_{T}^{part} (GeV/#it{c});#it{p}_{T}^{rec} (GeV/#it{c})", HistType::kTH2F, {axisParticlePt, axisParticlePt}); + addPairHistograms(directory + "/mc/particle/sameEvent", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + addPairHistograms(directory + "/mc/particle/mixedEvent", axisPairType, axisJetPt, axisDeltaEta, axisDeltaPhi, axisDeltaR); + addHistogram("/mc/particle/normalizationSameEvent", speciesName + " particle-level same-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); + addHistogram("/mc/particle/normalizationMixedEvent", speciesName + " particle-level mixed-event normalization;counter;weighted entries", HistType::kTH1F, {axisCounter}); + } + void addPairHistograms(std::string const& directory, AxisSpec const& axisPairType, AxisSpec const& axisJetPt, @@ -336,6 +411,18 @@ struct BaryonPairInJets { setNormalizationLabels(registry.get(HIST("normalization/mixedEventNominal"))); setNormalizationLabels(registry.get(HIST("normalization/mixedEventLoose"))); setNormalizationLabels(registry.get(HIST("normalization/mixedEventRandom"))); + if (output.fillPionControl.value) { + setNormalizationLabels(registry.get(HIST("control/pion/normalizationSameEvent"))); + setNormalizationLabels(registry.get(HIST("control/pion/normalizationMixedEvent"))); + setNormalizationLabels(registry.get(HIST("control/pion/mc/particle/normalizationSameEvent"))); + setNormalizationLabels(registry.get(HIST("control/pion/mc/particle/normalizationMixedEvent"))); + } + if (output.fillKaonControl.value) { + setNormalizationLabels(registry.get(HIST("control/kaon/normalizationSameEvent"))); + setNormalizationLabels(registry.get(HIST("control/kaon/normalizationMixedEvent"))); + setNormalizationLabels(registry.get(HIST("control/kaon/mc/particle/normalizationSameEvent"))); + setNormalizationLabels(registry.get(HIST("control/kaon/mc/particle/normalizationMixedEvent"))); + } } template @@ -356,7 +443,7 @@ struct BaryonPairInJets { if (fillCutFlow) { registry.fill(HIST("event/cutFlow"), 3.f); } - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, true, eventCuts.applyRctSelection.value)) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, true, eventCuts.applyRctSelection.value, static_cast(eventCuts.rctLabel))) { return false; } if (fillCutFlow) { @@ -436,10 +523,39 @@ struct BaryonPairInJets { return !isPionCompatible && !isKaonCompatible; } + bool passesControlPid(float pt, bool hasTof, float tpcNSigma, float tofNSigma, float tpcMaximum, float tofMaximum, float combinedMaximum) const + { + auto const pidMethod = static_cast(pidCuts.method.value); + tpcNSigma = std::abs(tpcNSigma); + tofNSigma = std::abs(tofNSigma); + if (pidMethod == ProtonPidMethod::TpcOnly || pt < pidCuts.tofPtThreshold.value) { + return tpcNSigma < tpcMaximum; + } + if (!hasTof) { + return false; + } + if (pidMethod == ProtonPidMethod::TpcTofCircular) { + return std::hypot(tpcNSigma, tofNSigma) < combinedMaximum; + } + return tpcNSigma < tpcMaximum && tofNSigma < tofMaximum; + } + + template + bool isSelectedPion(TrackType const& track) const + { + return passesControlPid(track.pt(), track.hasTOF(), track.tpcNSigmaPi(), track.tofNSigmaPi(), pidCuts.pionTpcNSigmaMax.value, pidCuts.pionTofNSigmaMax.value, pidCuts.pionCombinedNSigmaMax.value); + } + + template + bool isSelectedKaon(TrackType const& track) const + { + return passesControlPid(track.pt(), track.hasTOF(), track.tpcNSigmaKa(), track.tofNSigmaKa(), pidCuts.kaonTpcNSigmaMax.value, pidCuts.kaonTofNSigmaMax.value, pidCuts.kaonCombinedNSigmaMax.value); + } + template - std::vector collectBaryons(JetType const& jet, bool fillQa, float weight = 1.f) + SelectedParticles collectSelectedParticles(JetType const& jet, bool fillQa, float weight = 1.f, bool collectControls = true) { - std::vector baryons; + SelectedParticles selected; for (auto const& jetTrack : jet.template tracks_as()) { auto track = jetTrack.template track_as(); if (!isSelectedTrack(track)) { @@ -462,74 +578,250 @@ struct BaryonPairInJets { } } - float const rapidity = RecoDecay::y(std::array{track.px(), track.py(), track.pz()}, o2::constants::physics::MassProton); - if (std::abs(rapidity) > trackCuts.protonRapidityMax.value || !isSelectedProton(track)) { - continue; + SelectedParticle const particle{track.pt(), track.eta(), track.phi(), static_cast(track.sign())}; + if (isSelectedProton(track)) { + selected.baryons.push_back(particle); + if (fillQa) { + float const rapidity = RecoDecay::y(std::array{track.px(), track.py(), track.pz()}, o2::constants::physics::MassProton); + registry.fill(HIST("proton/pt"), particle.pt, weight); + registry.fill(HIST("proton/eta"), particle.eta, weight); + registry.fill(HIST("proton/rapidity"), rapidity, weight); + registry.fill(HIST("proton/tpcNSigmaVsPt"), track.pt(), track.tpcNSigmaPr(), weight); + if (track.hasTOF()) { + registry.fill(HIST("proton/tofNSigmaVsPt"), track.pt(), track.tofNSigmaPr(), weight); + } + registry.fill(HIST("proton/dcaXYVsPt"), track.pt(), track.dcaXY(), weight); + registry.fill(HIST("proton/dcaZVsPt"), track.pt(), track.dcaZ(), weight); + registry.fill(HIST("pid/selectedPtByCharge"), static_cast(particle.sign), particle.pt, weight); + } } - - baryons.push_back({track.pt(), track.eta(), track.phi(), rapidity, static_cast(track.sign())}); - if (fillQa) { - registry.fill(HIST("proton/pt"), track.pt(), weight); - registry.fill(HIST("proton/eta"), track.eta(), weight); - registry.fill(HIST("proton/rapidity"), rapidity, weight); - registry.fill(HIST("proton/tpcNSigmaVsPt"), track.pt(), track.tpcNSigmaPr(), weight); - if (track.hasTOF()) { - registry.fill(HIST("proton/tofNSigmaVsPt"), track.pt(), track.tofNSigmaPr(), weight); + if (collectControls && output.fillPionControl.value && isSelectedPion(track)) { + selected.pions.push_back(particle); + if (fillQa) { + registry.fill(HIST("control/pion/qa/pt"), particle.pt, weight); + registry.fill(HIST("control/pion/qa/tpcNSigmaVsPt"), particle.pt, track.tpcNSigmaPi(), weight); + if (track.hasTOF()) { + registry.fill(HIST("control/pion/qa/tofNSigmaVsPt"), particle.pt, track.tofNSigmaPi(), weight); + } + } + } + if (collectControls && output.fillKaonControl.value && isSelectedKaon(track)) { + selected.kaons.push_back(particle); + if (fillQa) { + registry.fill(HIST("control/kaon/qa/pt"), particle.pt, weight); + registry.fill(HIST("control/kaon/qa/tpcNSigmaVsPt"), particle.pt, track.tpcNSigmaKa(), weight); + if (track.hasTOF()) { + registry.fill(HIST("control/kaon/qa/tofNSigmaVsPt"), particle.pt, track.tofNSigmaKa(), weight); + } } - registry.fill(HIST("proton/dcaXYVsPt"), track.pt(), track.dcaXY(), weight); - registry.fill(HIST("proton/dcaZVsPt"), track.pt(), track.dcaZ(), weight); - registry.fill(HIST("pid/selectedPtByCharge"), static_cast(track.sign()), track.pt(), weight); } } - return baryons; + return selected; } - static PairType getPairType(SelectedBaryon const& baryon1, SelectedBaryon const& baryon2) + static PairType getPairType(SelectedParticle const& baryon1, SelectedParticle const& baryon2) { return baryon1.sign == baryon2.sign ? PairType::SameSign : PairType::OppositeSign; } - static ChargePairType getChargePairType(SelectedBaryon const& baryon1, SelectedBaryon const& baryon2) + static ChargePairType getChargePairType(SelectedParticle const& baryon1, SelectedParticle const& baryon2) { if (baryon1.sign > 0 && baryon2.sign > 0) { - return ChargePairType::ProtonProton; + return ChargePairType::PositivePositive; } if (baryon1.sign < 0 && baryon2.sign < 0) { - return ChargePairType::AntiprotonAntiproton; + return ChargePairType::NegativeNegative; + } + return ChargePairType::OppositeSign; + } + + template + void fillControlNormalization(bool isMixed, float counter, float weight) + { + if constexpr (Species == ControlSpecies::Pion) { + if (isMixed) { + registry.fill(HIST("control/pion/normalizationMixedEvent"), counter, weight); + } else { + registry.fill(HIST("control/pion/normalizationSameEvent"), counter, weight); + } + } else { + if (isMixed) { + registry.fill(HIST("control/kaon/normalizationMixedEvent"), counter, weight); + } else { + registry.fill(HIST("control/kaon/normalizationSameEvent"), counter, weight); + } + } + } + + template + void fillControlParticleNormalization(bool isMixed, float counter, float weight) + { + if constexpr (Species == ControlSpecies::Pion) { + if (isMixed) { + registry.fill(HIST("control/pion/mc/particle/normalizationMixedEvent"), counter, weight); + } else { + registry.fill(HIST("control/pion/mc/particle/normalizationSameEvent"), counter, weight); + } + } else { + if (isMixed) { + registry.fill(HIST("control/kaon/mc/particle/normalizationMixedEvent"), counter, weight); + } else { + registry.fill(HIST("control/kaon/mc/particle/normalizationSameEvent"), counter, weight); + } + } + } + + template + void fillControlPair(bool isMixed, PairType pairType, float jetPt, float deltaEta, float deltaPhi, float deltaR, float weight) + { + auto fillPairHistograms = [&](auto const& sameSignName, auto const& oppositeSignName, auto const& deltaRSameSignName, auto const& deltaROppositeSignName, auto const& sparseName) { + if (pairType == PairType::SameSign) { + registry.fill(sameSignName, deltaEta, deltaPhi, 0.5f * weight); + registry.fill(sameSignName, -deltaEta, -deltaPhi, 0.5f * weight); + registry.fill(deltaRSameSignName, deltaR, weight); + } else { + registry.fill(oppositeSignName, deltaEta, deltaPhi, 0.5f * weight); + registry.fill(oppositeSignName, -deltaEta, -deltaPhi, 0.5f * weight); + registry.fill(deltaROppositeSignName, deltaR, weight); + } + registry.fill(sparseName, static_cast(pairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(sparseName, static_cast(pairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + }; + + if constexpr (Species == ControlSpecies::Pion && !IsParticleLevel) { + if (isMixed) { + fillPairHistograms(HIST("control/pion/mixedEvent/sameSignDeltaEtaDeltaPhi"), HIST("control/pion/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), HIST("control/pion/mixedEvent/sameSignDeltaR"), HIST("control/pion/mixedEvent/oppositeSignDeltaR"), HIST("control/pion/mixedEvent/pairSparse")); + } else { + fillPairHistograms(HIST("control/pion/sameEvent/sameSignDeltaEtaDeltaPhi"), HIST("control/pion/sameEvent/oppositeSignDeltaEtaDeltaPhi"), HIST("control/pion/sameEvent/sameSignDeltaR"), HIST("control/pion/sameEvent/oppositeSignDeltaR"), HIST("control/pion/sameEvent/pairSparse")); + } + } else if constexpr (Species == ControlSpecies::Kaon && !IsParticleLevel) { + if (isMixed) { + fillPairHistograms(HIST("control/kaon/mixedEvent/sameSignDeltaEtaDeltaPhi"), HIST("control/kaon/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), HIST("control/kaon/mixedEvent/sameSignDeltaR"), HIST("control/kaon/mixedEvent/oppositeSignDeltaR"), HIST("control/kaon/mixedEvent/pairSparse")); + } else { + fillPairHistograms(HIST("control/kaon/sameEvent/sameSignDeltaEtaDeltaPhi"), HIST("control/kaon/sameEvent/oppositeSignDeltaEtaDeltaPhi"), HIST("control/kaon/sameEvent/sameSignDeltaR"), HIST("control/kaon/sameEvent/oppositeSignDeltaR"), HIST("control/kaon/sameEvent/pairSparse")); + } + } else if constexpr (Species == ControlSpecies::Pion) { + if (isMixed) { + fillPairHistograms(HIST("control/pion/mc/particle/mixedEvent/sameSignDeltaEtaDeltaPhi"), HIST("control/pion/mc/particle/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), HIST("control/pion/mc/particle/mixedEvent/sameSignDeltaR"), HIST("control/pion/mc/particle/mixedEvent/oppositeSignDeltaR"), HIST("control/pion/mc/particle/mixedEvent/pairSparse")); + } else { + fillPairHistograms(HIST("control/pion/mc/particle/sameEvent/sameSignDeltaEtaDeltaPhi"), HIST("control/pion/mc/particle/sameEvent/oppositeSignDeltaEtaDeltaPhi"), HIST("control/pion/mc/particle/sameEvent/sameSignDeltaR"), HIST("control/pion/mc/particle/sameEvent/oppositeSignDeltaR"), HIST("control/pion/mc/particle/sameEvent/pairSparse")); + } + } else { + if (isMixed) { + fillPairHistograms(HIST("control/kaon/mc/particle/mixedEvent/sameSignDeltaEtaDeltaPhi"), HIST("control/kaon/mc/particle/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), HIST("control/kaon/mc/particle/mixedEvent/sameSignDeltaR"), HIST("control/kaon/mc/particle/mixedEvent/oppositeSignDeltaR"), HIST("control/kaon/mc/particle/mixedEvent/pairSparse")); + } else { + fillPairHistograms(HIST("control/kaon/mc/particle/sameEvent/sameSignDeltaEtaDeltaPhi"), HIST("control/kaon/mc/particle/sameEvent/oppositeSignDeltaEtaDeltaPhi"), HIST("control/kaon/mc/particle/sameEvent/sameSignDeltaR"), HIST("control/kaon/mc/particle/sameEvent/oppositeSignDeltaR"), HIST("control/kaon/mc/particle/sameEvent/pairSparse")); + } + } + } + + template + void fillControlSameJetPairs(std::vector const& particles, float jetPt, float weight) + { + if constexpr (IsParticleLevel) { + fillControlParticleNormalization(false, 2.f, weight); + } else { + fillControlNormalization(false, 2.f, weight); + } + for (size_t first = 0; first < particles.size(); ++first) { + for (size_t second = first + 1; second < particles.size(); ++second) { + PairType const pairType = getPairType(particles[first], particles[second]); + float const deltaEta = particles[first].eta - particles[second].eta; + float const deltaPhi = RecoDecay::constrainAngle(particles[first].phi - particles[second].phi, -o2::constants::math::PI); + float const deltaR = std::hypot(deltaEta, deltaPhi); + fillControlPair(false, pairType, jetPt, deltaEta, deltaPhi, deltaR, weight); + if constexpr (IsParticleLevel) { + fillControlParticleNormalization(false, pairType == PairType::SameSign ? 3.f : 4.f, weight); + } else { + fillControlNormalization(false, pairType == PairType::SameSign ? 3.f : 4.f, weight); + } + } + } + } + + template + void fillControlMixedJetPairs(std::vector const& particles1, + std::vector const& particles2, + JetType1 const& jet1, + JetType2 const& jet2, + float weight) + { + if constexpr (IsParticleLevel) { + fillControlParticleNormalization(true, 2.f, weight); + } else { + fillControlNormalization(true, 2.f, weight); + } + float const mixedJetPt = 0.5f * (jet1.pt() + jet2.pt()); + for (auto const& particle1 : particles1) { + float const localEta1 = particle1.eta - jet1.eta(); + float const localPhi1 = RecoDecay::constrainAngle(particle1.phi - jet1.phi(), -o2::constants::math::PI); + for (auto const& particle2 : particles2) { + float const localEta2 = particle2.eta - jet2.eta(); + float const localPhi2 = RecoDecay::constrainAngle(particle2.phi - jet2.phi(), -o2::constants::math::PI); + float const deltaEta = localEta1 - localEta2; + float const deltaPhi = RecoDecay::constrainAngle(localPhi1 - localPhi2, -o2::constants::math::PI); + float const deltaR = std::hypot(deltaEta, deltaPhi); + PairType const pairType = getPairType(particle1, particle2); + fillControlPair(true, pairType, mixedJetPt, deltaEta, deltaPhi, deltaR, weight); + if constexpr (IsParticleLevel) { + fillControlParticleNormalization(true, pairType == PairType::SameSign ? 3.f : 4.f, weight); + } else { + fillControlNormalization(true, pairType == PairType::SameSign ? 3.f : 4.f, weight); + } + } } - return ChargePairType::ProtonAntiproton; } template - std::vector collectInclusiveBaryons(TracksType const& tracks) + std::vector collectInclusiveBaryons(TracksType const& tracks) { - std::vector baryons; + std::vector baryons; for (auto const& track : tracks) { if (!isSelectedTrack(track)) { continue; } - float const rapidity = RecoDecay::y(std::array{track.px(), track.py(), track.pz()}, o2::constants::physics::MassProton); - if (std::abs(rapidity) > trackCuts.protonRapidityMax.value || !isSelectedProton(track)) { + if (!isSelectedProton(track)) { continue; } - baryons.push_back({track.pt(), track.eta(), track.phi(), rapidity, static_cast(track.sign())}); + baryons.push_back({track.pt(), track.eta(), track.phi(), static_cast(track.sign())}); } return baryons; } template - bool isSelectedInclusiveCollision(CollisionType const& collision) const + bool isSelectedInclusiveCollision(CollisionType const& collision) { - return std::abs(collision.posZ()) <= eventCuts.vertexZMax.value && collision.sel8() && - collision.selection_bit(aod::evsel::kNoITSROFrameBorder) && - collision.selection_bit(aod::evsel::kNoTimeFrameBorder) && - collision.selection_bit(aod::evsel::kNoSameBunchPileup) && - collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV); + if (std::abs(collision.posZ()) > eventCuts.vertexZMax.value || + collision.trackOccupancyInTimeRange() < eventCuts.occupancyMin.value || + collision.trackOccupancyInTimeRange() > eventCuts.occupancyMax.value) { + return false; + } + + uint16_t const collisionSelection = jetderiveddatautilities::setEventSelectionBit(collision); + bool isOrCondition = false; + bool passedOrCondition = false; + for (auto const& selectionBit : eventSelectionBits) { + if (selectionBit == -1) { + isOrCondition = true; + continue; + } + bool const passed = collisionSelection & (1U << selectionBit); + if (!isOrCondition && !passed) { + return false; + } + if (isOrCondition && passed) { + passedOrCondition = true; + } + } + if (isOrCondition && !passedOrCondition) { + return false; + } + return !eventCuts.applyRctSelection.value || inclusiveRctChecker.checkTable(collision); } void fillInclusivePair(bool isMixed, - SelectedBaryon const& baryon1, - SelectedBaryon const& baryon2, + SelectedParticle const& baryon1, + SelectedParticle const& baryon2, float weight) { PairType const pairType = getPairType(baryon1, baryon2); @@ -545,8 +837,10 @@ struct BaryonPairInJets { registry.fill(HIST("inclusive/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); registry.fill(HIST("inclusive/mixedEvent/oppositeSignDeltaEtaDeltaPhi"), -deltaEta, reflectedDeltaPhi, 0.5f * weight); } - registry.fill(HIST("inclusive/mixedEvent/chargePairSparse"), static_cast(chargePairType), deltaEta, deltaPhi, 0.5f * weight); - registry.fill(HIST("inclusive/mixedEvent/chargePairSparse"), static_cast(chargePairType), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + if (output.fillChargeSeparated.value) { + registry.fill(HIST("inclusive/mixedEvent/chargePairSparse"), static_cast(chargePairType), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("inclusive/mixedEvent/chargePairSparse"), static_cast(chargePairType), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + } } else { if (pairType == PairType::SameSign) { registry.fill(HIST("inclusive/sameEvent/sameSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); @@ -555,8 +849,10 @@ struct BaryonPairInJets { registry.fill(HIST("inclusive/sameEvent/oppositeSignDeltaEtaDeltaPhi"), deltaEta, deltaPhi, 0.5f * weight); registry.fill(HIST("inclusive/sameEvent/oppositeSignDeltaEtaDeltaPhi"), -deltaEta, reflectedDeltaPhi, 0.5f * weight); } - registry.fill(HIST("inclusive/sameEvent/chargePairSparse"), static_cast(chargePairType), deltaEta, deltaPhi, 0.5f * weight); - registry.fill(HIST("inclusive/sameEvent/chargePairSparse"), static_cast(chargePairType), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + if (output.fillChargeSeparated.value) { + registry.fill(HIST("inclusive/sameEvent/chargePairSparse"), static_cast(chargePairType), deltaEta, deltaPhi, 0.5f * weight); + registry.fill(HIST("inclusive/sameEvent/chargePairSparse"), static_cast(chargePairType), -deltaEta, reflectedDeltaPhi, 0.5f * weight); + } } } @@ -599,8 +895,10 @@ struct BaryonPairInJets { } else { registry.fill(HIST("mixedEvent/nominal/oppositeSignDeltaR"), deltaR, weight); } - registry.fill(HIST("mixedEvent/nominal/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); - registry.fill(HIST("mixedEvent/nominal/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + if (output.fillChargeSeparated.value) { + registry.fill(HIST("mixedEvent/nominal/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("mixedEvent/nominal/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } } else if constexpr (Reference == MixedReference::Loose) { fillSymmetrized(HIST("mixedEvent/loose/sameSignDeltaEtaDeltaPhi"), HIST("mixedEvent/loose/oppositeSignDeltaEtaDeltaPhi"), HIST("mixedEvent/loose/pairSparse")); if (pairType == PairType::SameSign) { @@ -608,8 +906,10 @@ struct BaryonPairInJets { } else { registry.fill(HIST("mixedEvent/loose/oppositeSignDeltaR"), deltaR, weight); } - registry.fill(HIST("mixedEvent/loose/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); - registry.fill(HIST("mixedEvent/loose/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + if (output.fillChargeSeparated.value) { + registry.fill(HIST("mixedEvent/loose/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("mixedEvent/loose/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } } else { fillSymmetrized(HIST("mixedEvent/random/sameSignDeltaEtaDeltaPhi"), HIST("mixedEvent/random/oppositeSignDeltaEtaDeltaPhi"), HIST("mixedEvent/random/pairSparse")); if (pairType == PairType::SameSign) { @@ -617,8 +917,10 @@ struct BaryonPairInJets { } else { registry.fill(HIST("mixedEvent/random/oppositeSignDeltaR"), deltaR, weight); } - registry.fill(HIST("mixedEvent/random/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); - registry.fill(HIST("mixedEvent/random/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + if (output.fillChargeSeparated.value) { + registry.fill(HIST("mixedEvent/random/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("mixedEvent/random/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } } } @@ -648,29 +950,45 @@ struct BaryonPairInJets { } registry.fill(HIST("sameEvent/pairSparse"), pairTypeValue, jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); registry.fill(HIST("sameEvent/pairSparse"), pairTypeValue, jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); - registry.fill(HIST("sameEvent/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); - registry.fill(HIST("sameEvent/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + if (output.fillChargeSeparated.value) { + registry.fill(HIST("sameEvent/chargePairSparse"), static_cast(chargePairType), jetPt, deltaEta, deltaPhi, deltaR, 0.5f * weight); + registry.fill(HIST("sameEvent/chargePairSparse"), static_cast(chargePairType), jetPt, -deltaEta, -deltaPhi, deltaR, 0.5f * weight); + } } template - std::vector collectParticleBaryons(JetType const& jet, bool fillQa, float weight) + SelectedParticles collectParticleCandidates(JetType const& jet, bool fillQa, float weight) { - std::vector baryons; + SelectedParticles selected; for (auto const& particle : jet.template tracks_as()) { - if (std::abs(particle.pdgCode()) != PDG_t::kProton || !particle.isPhysicalPrimary()) { + if (!particle.isPhysicalPrimary()) { continue; } if (particle.pt() < trackCuts.protonPtMin.value || particle.pt() >= trackCuts.protonPtMax.value || - std::abs(particle.eta()) > trackCuts.protonEtaMax.value || std::abs(particle.y()) > trackCuts.protonRapidityMax.value) { + std::abs(particle.eta()) > trackCuts.protonEtaMax.value) { continue; } int8_t const sign = particle.pdgCode() > 0 ? 1 : -1; - baryons.push_back({particle.pt(), particle.eta(), particle.phi(), particle.y(), sign}); - if (fillQa) { - registry.fill(HIST("mc/pid/generatedPrimaryProtonPt"), particle.pt(), weight); + SelectedParticle const selectedParticle{particle.pt(), particle.eta(), particle.phi(), sign}; + int const absolutePdg = std::abs(particle.pdgCode()); + if (absolutePdg == PDG_t::kProton) { + selected.baryons.push_back(selectedParticle); + if (fillQa) { + registry.fill(HIST("mc/pid/generatedPrimaryProtonPt"), selectedParticle.pt, weight); + } + } else if (output.fillPionControl.value && absolutePdg == PDG_t::kPiPlus) { + selected.pions.push_back(selectedParticle); + if (fillQa) { + registry.fill(HIST("control/pion/mc/generatedPrimaryPt"), selectedParticle.pt, weight); + } + } else if (output.fillKaonControl.value && absolutePdg == PDG_t::kKPlus) { + selected.kaons.push_back(selectedParticle); + if (fillQa) { + registry.fill(HIST("control/kaon/mc/generatedPrimaryPt"), selectedParticle.pt, weight); + } } } - return baryons; + return selected; } void fillParticlePair(bool isMixed, PairType pairType, float jetPt, float deltaEta, float deltaPhi, float deltaR, float weight) @@ -702,6 +1020,48 @@ struct BaryonPairInJets { } } + template + void fillControlMCPidQa(TrackType const& track, ParticleType const& particle, float weight) + { + int const expectedPdg = Species == ControlSpecies::Pion ? PDG_t::kPiPlus : PDG_t::kKPlus; + bool const isTruePrimary = std::abs(particle.pdgCode()) == expectedPdg && particle.isPhysicalPrimary(); + bool const isSelected = Species == ControlSpecies::Pion ? isSelectedPion(track) : isSelectedKaon(track); + int origin = 2; + if (particle.isPhysicalPrimary()) { + origin = 1; + } else if (particle.getProcess() == TMCProcess::kPDecay) { + origin = 3; + } + + if constexpr (Species == ControlSpecies::Pion) { + if (isTruePrimary) { + registry.fill(HIST("control/pion/mc/reconstructedTruePt"), track.pt(), weight); + } + if (!isSelected) { + return; + } + registry.fill(HIST("control/pion/mc/selectedCandidateOriginVsPt"), track.pt(), static_cast(origin), weight); + registry.fill(HIST("control/pion/mc/selectedCandidateTruePdgVsPt"), track.pt(), static_cast(particle.pdgCode()), weight); + if (isTruePrimary) { + registry.fill(HIST("control/pion/mc/selectedTruePt"), track.pt(), weight); + registry.fill(HIST("control/pion/mc/ptResponsePrimary"), particle.pt(), track.pt(), weight); + } + } else { + if (isTruePrimary) { + registry.fill(HIST("control/kaon/mc/reconstructedTruePt"), track.pt(), weight); + } + if (!isSelected) { + return; + } + registry.fill(HIST("control/kaon/mc/selectedCandidateOriginVsPt"), track.pt(), static_cast(origin), weight); + registry.fill(HIST("control/kaon/mc/selectedCandidateTruePdgVsPt"), track.pt(), static_cast(particle.pdgCode()), weight); + if (isTruePrimary) { + registry.fill(HIST("control/kaon/mc/selectedTruePt"), track.pt(), weight); + registry.fill(HIST("control/kaon/mc/ptResponsePrimary"), particle.pt(), track.pt(), weight); + } + } + } + template void fillMCPidQa(JetType const& jet, float weight) { @@ -711,12 +1071,17 @@ struct BaryonPairInJets { continue; } auto particle = track.template mcParticle_as(); + if (output.fillPionControl.value) { + fillControlMCPidQa(track, particle, weight); + } + if (output.fillKaonControl.value) { + fillControlMCPidQa(track, particle, weight); + } bool const isTruePrimaryProton = std::abs(particle.pdgCode()) == PDG_t::kProton && particle.isPhysicalPrimary(); if (isTruePrimaryProton) { registry.fill(HIST("mc/pid/reconstructedTrueProtonPt"), track.pt(), weight); } - float const rapidity = RecoDecay::y(std::array{track.px(), track.py(), track.pz()}, o2::constants::physics::MassProton); - if (std::abs(rapidity) > trackCuts.protonRapidityMax.value || !isSelectedProton(track)) { + if (!isSelectedProton(track)) { continue; } int origin = 2; @@ -741,12 +1106,20 @@ struct BaryonPairInJets { return; } - registry.fill(HIST("event/zVertex"), collision.posZ(), weight); - registry.fill(HIST("event/multiplicityFT0C"), collision.multFT0C(), weight); + if (output.fillQa.value) { + registry.fill(HIST("event/zVertex"), collision.posZ(), weight); + registry.fill(HIST("event/multiplicityFT0C"), collision.multFT0C(), weight); + } registry.fill(HIST("normalization/sameEvent"), 1.f, weight); + if (output.fillPionControl.value) { + fillControlNormalization(false, 1.f, weight); + } + if (output.fillKaonControl.value) { + fillControlNormalization(false, 1.f, weight); + } for (auto const& jet : jets) { - if (isSelectedJetGeometry(jet)) { + if (output.fillQa.value && isSelectedJetGeometry(jet)) { registry.fill(HIST("jet/ptBeforePtSelection"), jet.pt(), weight); } if (!isSelectedJet(jet)) { @@ -754,23 +1127,28 @@ struct BaryonPairInJets { } registry.fill(HIST("normalization/sameEvent"), 2.f, weight); - registry.fill(HIST("jet/pt"), jet.pt(), weight); - registry.fill(HIST("jet/eta"), jet.eta(), weight); - registry.fill(HIST("jet/phi"), jet.phi(), weight); - registry.fill(HIST("jet/nConstituents"), static_cast(jet.tracksIds().size()), weight); - registry.fill(HIST("jet/nConstituentsVsPt"), jet.pt(), static_cast(jet.tracksIds().size()), weight); - - auto baryons = collectBaryons(jet, true, weight); - int numberOfProtons = 0; - int numberOfAntiprotons = 0; - for (auto const& baryon : baryons) { - baryon.sign > 0 ? ++numberOfProtons : ++numberOfAntiprotons; - } - registry.fill(HIST("jet/nSelectedBaryons"), static_cast(baryons.size()), weight); - registry.fill(HIST("jet/nProtons"), numberOfProtons, weight); - registry.fill(HIST("jet/nAntiprotons"), numberOfAntiprotons, weight); - registry.fill(HIST("jet/nProtonsVsAntiprotons"), numberOfProtons, numberOfAntiprotons, weight); - registry.fill(HIST("jet/nBaryonsVsPt"), jet.pt(), static_cast(baryons.size()), weight); + if (output.fillQa.value) { + registry.fill(HIST("jet/pt"), jet.pt(), weight); + registry.fill(HIST("jet/eta"), jet.eta(), weight); + registry.fill(HIST("jet/phi"), jet.phi(), weight); + registry.fill(HIST("jet/nConstituents"), static_cast(jet.tracksIds().size()), weight); + registry.fill(HIST("jet/nConstituentsVsPt"), jet.pt(), static_cast(jet.tracksIds().size()), weight); + } + + auto selectedParticles = collectSelectedParticles(jet, output.fillQa.value, weight); + auto const& baryons = selectedParticles.baryons; + if (output.fillQa.value) { + int numberOfProtons = 0; + int numberOfAntiprotons = 0; + for (auto const& baryon : baryons) { + baryon.sign > 0 ? ++numberOfProtons : ++numberOfAntiprotons; + } + registry.fill(HIST("jet/nSelectedBaryons"), static_cast(baryons.size()), weight); + registry.fill(HIST("jet/nProtons"), numberOfProtons, weight); + registry.fill(HIST("jet/nAntiprotons"), numberOfAntiprotons, weight); + registry.fill(HIST("jet/nProtonsVsAntiprotons"), numberOfProtons, numberOfAntiprotons, weight); + registry.fill(HIST("jet/nBaryonsVsPt"), jet.pt(), static_cast(baryons.size()), weight); + } for (size_t first = 0; first < baryons.size(); ++first) { for (size_t second = first + 1; second < baryons.size(); ++second) { @@ -780,10 +1158,12 @@ struct BaryonPairInJets { float const deltaPhi = RecoDecay::constrainAngle(baryons[first].phi - baryons[second].phi, -o2::constants::math::PI); float const deltaR = std::hypot(deltaEta, deltaPhi); fillSameEventPair(pairType, chargePairType, jet.pt(), deltaEta, deltaPhi, deltaR, weight); - if (pairType == PairType::SameSign) { - registry.fill(HIST("jet/sameSignPairsVsPt"), jet.pt(), weight); - } else { - registry.fill(HIST("jet/oppositeSignPairsVsPt"), jet.pt(), weight); + if (output.fillQa.value) { + if (pairType == PairType::SameSign) { + registry.fill(HIST("jet/sameSignPairsVsPt"), jet.pt(), weight); + } else { + registry.fill(HIST("jet/oppositeSignPairsVsPt"), jet.pt(), weight); + } } float const counterBin = pairType == PairType::SameSign ? 3.f : 4.f; @@ -793,6 +1173,12 @@ struct BaryonPairInJets { } } } + if (output.fillPionControl.value) { + fillControlSameJetPairs(selectedParticles.pions, jet.pt(), weight); + } + if (output.fillKaonControl.value) { + fillControlSameJetPairs(selectedParticles.kaons, jet.pt(), weight); + } } } @@ -821,8 +1207,11 @@ struct BaryonPairInJets { registry.fill(HIST("mixedEvent/jetNchDifference"), jetNchDifference, weight); } - auto baryons1 = collectBaryons(jet1, false); - auto baryons2 = collectBaryons(jet2, false); + constexpr bool CollectControls = Reference == MixedReference::Nominal; + auto selectedParticles1 = collectSelectedParticles(jet1, false, 1.f, CollectControls); + auto selectedParticles2 = collectSelectedParticles(jet2, false, 1.f, CollectControls); + auto const& baryons1 = selectedParticles1.baryons; + auto const& baryons2 = selectedParticles2.baryons; float const mixedJetPt = 0.5f * (jet1.pt() + jet2.pt()); fillMixedNormalization(2.f, weight); @@ -844,6 +1233,14 @@ struct BaryonPairInJets { } } } + if constexpr (Reference == MixedReference::Nominal) { + if (output.fillPionControl.value) { + fillControlMixedJetPairs(selectedParticles1.pions, selectedParticles2.pions, jet1, jet2, weight); + } + if (output.fillKaonControl.value) { + fillControlMixedJetPairs(selectedParticles1.kaons, selectedParticles2.kaons, jet1, jet2, weight); + } + } } void processMixedEvent(aod::JetCollisions const& collisions, @@ -864,8 +1261,16 @@ struct BaryonPairInJets { float const weight = 1.f / static_cast(iterator.currentWindowNeighbours()); registry.fill(HIST("normalization/mixedEventNominal"), 1.f, weight); - registry.fill(HIST("normalization/mixedEventLoose"), 1.f, weight); - registry.fill(HIST("normalization/mixedEventRandom"), 1.f, weight); + if (output.fillPionControl.value) { + fillControlNormalization(true, 1.f, weight); + } + if (output.fillKaonControl.value) { + fillControlNormalization(true, 1.f, weight); + } + if (output.fillReferenceVariations.value) { + registry.fill(HIST("normalization/mixedEventLoose"), 1.f, weight); + registry.fill(HIST("normalization/mixedEventRandom"), 1.f, weight); + } for (auto const& jet1 : jets1) { if (!isSelectedJet(jet1)) { @@ -876,8 +1281,10 @@ struct BaryonPairInJets { continue; } fillMixedJetPair(jet1, jet2, weight); - fillMixedJetPair(jet1, jet2, weight); - fillMixedJetPair(jet1, jet2, weight); + if (output.fillReferenceVariations.value) { + fillMixedJetPair(jet1, jet2, weight); + fillMixedJetPair(jet1, jet2, weight); + } } } } @@ -923,8 +1330,16 @@ struct BaryonPairInJets { } float const weight = collision1.mcCollision().weight() / static_cast(iterator.currentWindowNeighbours()); registry.fill(HIST("normalization/mixedEventNominal"), 1.f, weight); - registry.fill(HIST("normalization/mixedEventLoose"), 1.f, weight); - registry.fill(HIST("normalization/mixedEventRandom"), 1.f, weight); + if (output.fillPionControl.value) { + fillControlNormalization(true, 1.f, weight); + } + if (output.fillKaonControl.value) { + fillControlNormalization(true, 1.f, weight); + } + if (output.fillReferenceVariations.value) { + registry.fill(HIST("normalization/mixedEventLoose"), 1.f, weight); + registry.fill(HIST("normalization/mixedEventRandom"), 1.f, weight); + } for (auto const& jet1 : jets1) { if (!isSelectedJet(jet1)) { continue; @@ -934,8 +1349,10 @@ struct BaryonPairInJets { continue; } fillMixedJetPair(jet1, jet2, weight); - fillMixedJetPair(jet1, jet2, weight); - fillMixedJetPair(jet1, jet2, weight); + if (output.fillReferenceVariations.value) { + fillMixedJetPair(jet1, jet2, weight); + fillMixedJetPair(jet1, jet2, weight); + } } } } @@ -951,13 +1368,20 @@ struct BaryonPairInJets { } float const weight = collision.weight(); registry.fill(HIST("mc/particle/normalizationSameEvent"), 1.f, weight); + if (output.fillPionControl.value) { + fillControlParticleNormalization(false, 1.f, weight); + } + if (output.fillKaonControl.value) { + fillControlParticleNormalization(false, 1.f, weight); + } for (auto const& jet : jets) { if (!isSelectedJet(jet)) { continue; } registry.fill(HIST("mc/jet/allParticlePt"), jet.pt(), weight); registry.fill(HIST("mc/particle/normalizationSameEvent"), 2.f, weight); - auto baryons = collectParticleBaryons(jet, true, weight); + auto selectedParticles = collectParticleCandidates(jet, true, weight); + auto const& baryons = selectedParticles.baryons; for (size_t first = 0; first < baryons.size(); ++first) { for (size_t second = first + 1; second < baryons.size(); ++second) { PairType const pairType = getPairType(baryons[first], baryons[second]); @@ -968,6 +1392,12 @@ struct BaryonPairInJets { registry.fill(HIST("mc/particle/normalizationSameEvent"), pairType == PairType::SameSign ? 3.f : 4.f, weight); } } + if (output.fillPionControl.value) { + fillControlSameJetPairs(selectedParticles.pions, jet.pt(), weight); + } + if (output.fillKaonControl.value) { + fillControlSameJetPairs(selectedParticles.kaons, jet.pt(), weight); + } } } PROCESS_SWITCH(BaryonPairInJets, processMCP, "Particle-level MC baryon pairs", false); @@ -987,6 +1417,12 @@ struct BaryonPairInJets { } float const weight = collision1.weight() / static_cast(iterator.currentWindowNeighbours()); registry.fill(HIST("mc/particle/normalizationMixedEvent"), 1.f, weight); + if (output.fillPionControl.value) { + fillControlParticleNormalization(true, 1.f, weight); + } + if (output.fillKaonControl.value) { + fillControlParticleNormalization(true, 1.f, weight); + } for (auto const& jet1 : jets1) { if (!isSelectedJet(jet1)) { continue; @@ -995,8 +1431,10 @@ struct BaryonPairInJets { if (!isSelectedJet(jet2)) { continue; } - auto baryons1 = collectParticleBaryons(jet1, false, weight); - auto baryons2 = collectParticleBaryons(jet2, false, weight); + auto selectedParticles1 = collectParticleCandidates(jet1, false, weight); + auto selectedParticles2 = collectParticleCandidates(jet2, false, weight); + auto const& baryons1 = selectedParticles1.baryons; + auto const& baryons2 = selectedParticles2.baryons; registry.fill(HIST("mc/particle/normalizationMixedEvent"), 2.f, weight); for (auto const& baryon1 : baryons1) { float const localEta1 = baryon1.eta - jet1.eta(); @@ -1012,6 +1450,12 @@ struct BaryonPairInJets { registry.fill(HIST("mc/particle/normalizationMixedEvent"), pairType == PairType::SameSign ? 3.f : 4.f, weight); } } + if (output.fillPionControl.value) { + fillControlMixedJetPairs(selectedParticles1.pions, selectedParticles2.pions, jet1, jet2, weight); + } + if (output.fillKaonControl.value) { + fillControlMixedJetPairs(selectedParticles1.kaons, selectedParticles2.kaons, jet1, jet2, weight); + } } } } @@ -1063,7 +1507,7 @@ struct BaryonPairInJets { if (!isSelectedInclusiveCollision(collision)) { return; } - registry.fill(HIST("inclusive/event/multiplicity"), collision.multNTracksPV()); + registry.fill(HIST("inclusive/event/multiplicity"), collision.multFT0C()); registry.fill(HIST("inclusive/normalization/sameEvent"), 1.f); auto baryons = collectInclusiveBaryons(tracks); registry.fill(HIST("inclusive/normalization/sameEvent"), 2.f, static_cast(baryons.size()));