diff --git a/PWGJE/Tasks/jetShape.cxx b/PWGJE/Tasks/jetShape.cxx index 0df2ae24253..c18834ab478 100644 --- a/PWGJE/Tasks/jetShape.cxx +++ b/PWGJE/Tasks/jetShape.cxx @@ -24,6 +24,7 @@ #include "Common/DataModel/TrackSelectionTables.h" #include +#include #include #include #include @@ -36,6 +37,7 @@ #include #include +#include #include #include #include @@ -47,6 +49,10 @@ using namespace o2::framework::expressions; struct JetShapeTask { + static constexpr float PionMass = static_cast(o2::constants::physics::MassPionCharged); + static constexpr float ProtonMass = static_cast(o2::constants::physics::MassProton); + static constexpr float RandomConeNegativeSideProbability = 0.5f; + Configurable nBinsNSigma{"nBinsNSigma", 100, "Number of nsigma bins"}; Configurable nSigmaMin{"nSigmaMin", -10.0f, "Min value of nsigma"}; Configurable nSigmaMax{"nSigmaMax", 10.0f, "Max value of nsigma"}; @@ -75,12 +81,19 @@ struct JetShapeTask { Configurable nBinsDistance{"nBinsDistance", 7, "Number of distance bins"}; Configurable distanceMax{"distanceMax", 0.7f, "Max value of distance"}; Configurable nSigmaTofCut{"nSigmaTofCut", 2.0f, "Number of sigma cut for TOF PID"}; - Configurable tpcNSigmaPrMin{"tpcNSigmaPrMin", -3.5f, "Min value of tpcNsigmaProton"}; - Configurable tpcNSigmaPrMax{"tpcNSigmaPrMax", 0.5f, "Max value of tpcNsigmaProton"}; - Configurable tpcNSigmaPiMin{"tpcNSigmaPiMin", -0.5f, "Min value of tpcNsigmaPion"}; - Configurable tpcNSigmaPiMax{"tpcNSigmaPiMax", 3.5f, "Max value of tpcNsigmaPion"}; - Configurable randomConeDeltaPhiMin{"randomConeDeltaPhiMin", static_cast(o2::constants::math::PIThird), "Minimum delta phi for random cone"}; - Configurable randomConeDeltaPhiMax{"randomConeDeltaPhiMax", static_cast(2.0f * o2::constants::math::PIThird), "Maximum delta phi for random cone"}; + Configurable tpcNSigmaPrMin{"tpcNSigmaPrMin", -2.0f, "Min value of tpcNsigmaProton"}; + Configurable tpcNSigmaPrMax{"tpcNSigmaPrMax", 2.0f, "Max value of tpcNsigmaProton"}; + Configurable tpcNSigmaPiMin{"tpcNSigmaPiMin", -2.0f, "Min value of tpcNsigmaPion"}; + Configurable tpcNSigmaPiMax{"tpcNSigmaPiMax", 2.0f, "Max value of tpcNsigmaPion"}; + Configurable randomConeDeltaPhiMin{"randomConeDeltaPhiMin", static_cast(o2::constants::math::PIThird), "Minimum delta phi used only when randomConeFullAzimuth=false"}; + Configurable randomConeDeltaPhiMax{"randomConeDeltaPhiMax", static_cast(2.0f * o2::constants::math::PIThird), "Maximum delta phi used only when randomConeFullAzimuth=false"}; + Configurable enableRandomCone{"enableRandomCone", true, "Fill one accepted random cone per selected jet"}; + Configurable randomConeFullAzimuth{"randomConeFullAzimuth", true, "Generate random cones over the full azimuth instead of the old perpendicular-side window"}; + Configurable randomConeMaxAttempts{"randomConeMaxAttempts", 100, "Maximum attempts for each random cone"}; + Configurable randomConeSeed{"randomConeSeed", 0u, "TRandom3 seed; 0 lets ROOT choose a unique seed"}; + Configurable fillDetailedQA{"fillDetailedQA", false, "Fill high-volume dE/dx, beta and distance QA histograms"}; + Configurable applyCentralityCut{"applyCentralityCut", true, "Apply centralityMinForCut <= centrality < centralityMaxForCut inside this task"}; + Configurable useRapiditySelection{"useRapiditySelection", false, "Use species-dependent |y| cut in addition to the detector |eta| acceptance (cross-check only)"}; Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; @@ -89,7 +102,8 @@ struct JetShapeTask { Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; - Configurable trackSelections{"trackSelections", "globalTracks", "set track selections"}; + Configurable trackSelections{"trackSelections", "globalTracks", "nominal track selection for yields and tracking efficiency"}; + Configurable trackSelectionsForDca{"trackSelectionsForDca", "QualityTracksWDCA", "wide-DCA track selection used only for data DCA distributions and MC DCA templates"}; Configurable jetAreaFractionMin{"jetAreaFractionMin", -99.0, "used to make a cut on the jet areas"}; Configurable leadingConstituentPtMin{"leadingConstituentPtMin", 5.0, "minimum pT selection on jet constituent"}; @@ -100,13 +114,14 @@ struct JetShapeTask { // for ppi production Configurable etaTrUp{"etaTrUp", 0.9f, "maximum track eta"}; + Configurable trackPtMinForAnalysis{"trackPtMinForAnalysis", 0.15f, "minimum reconstructed and generated particle pT used for yields and efficiencies"}; Configurable dcaxyCutMax{"dcaxyCutMax", 2.0f, "maximum DCA xy"}; Configurable chi2ItsMax{"chi2ItsMax", 15.0f, "its chi2 cut"}; Configurable chi2TpcMax{"chi2TpcMax", 4.0f, "tpc chi2 cut"}; Configurable nclItsMin{"nclItsMin", 2.0f, "its # of cluster cut"}; Configurable nclTpcMin{"nclTpcMin", 70.0f, "tpc # if cluster cut"}; Configurable nclcrossTpcMin{"nclcrossTpcMin", 70.0f, "tpc # of crossedRows cut"}; - Configurable mcRapidityMax{"mcRapidityMax", 0.5f, "maximum mctrack y"}; + Configurable mcRapidityMax{"mcRapidityMax", 0.5f, "maximum |y| used only when useRapiditySelection=true"}; Configurable epsilon{"epsilon", 1e-6, "standard for aboid division of zero"}; Configurable maxDeltaEtaSafe{"maxDeltaEtaSafe", 0.9f, "maximum track eta for cut"}; Configurable nSigmaMaxForDcaxy{"nSigmaMaxForDcaxy", 4.0f, "maximum nSigma for DCAxy"}; @@ -114,16 +129,20 @@ struct JetShapeTask { Configurable triggerMasks{"triggerMasks", "", "possible JE Trigger masks: fJetChLowPt,fJetChHighPt,fTrackLowPt,fTrackHighPt,fJetD0ChLowPt,fJetD0ChHighPt,fJetLcChLowPt,fJetLcChHighPt,fEMCALReadout,fJetFullHighPt,fJetFullLowPt,fJetNeutralHighPt,fJetNeutralLowPt,fGammaVeryHighPtEMCAL,fGammaVeryHighPtDCAL,fGammaHighPtEMCAL,fGammaHighPtDCAL,fGammaLowPtEMCAL,fGammaLowPtDCAL,fGammaVeryLowPtEMCAL,fGammaVeryLowPtDCAL"}; HistogramRegistry registry{"registry"}; + TRandom3 randomNumber; std::vector eventSelectionBits; int trackSelection = -1; + int trackSelectionForDca = -1; std::vector triggerMaskBits; void init(o2::framework::InitContext&) { eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); + trackSelectionForDca = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelectionsForDca)); triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(triggerMasks); + randomNumber.SetSeed(randomConeSeed); // Histograms definition registry.add("tpcTofPi", "tpcTofPi", HistType::kTHnSparseD, {{nBinsPForCut.value, 0, pMax.value}, {nBinsNSigma.value, nSigmaMin.value, nSigmaMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); @@ -136,6 +155,8 @@ struct JetShapeTask { registry.add("tofBeta", "tofBeta", HistType::kTHnSparseD, {{nBinsPForBeta.value, 0, pMax.value}, {nBinsTofBeta.value, 0.4, 1.1}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("pVsPtForPr", "pVsPtForPr", HistType::kTHnSparseD, {{nBinsP.value, 0, pMax.value}, {nBinsPt.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("pVsPtForPi", "pVsPtPi", HistType::kTHnSparseD, {{nBinsP.value, 0, pMax.value}, {nBinsPt.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("dcaPiInc", "inclusive pion-candidate DCAxy from wide-DCA tracks", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsDcaxyForData.value, dcaxyMin.value, dcaxyMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("dcaPrInc", "inclusive proton-candidate DCAxy from wide-DCA tracks", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsDcaxyForData.value, dcaxyMin.value, dcaxyMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("trackPhi", "trackPhi", HistType::kTH1F, {{80, -1, 7}}); registry.add("trackEta", "trackEta", HistType::kTH1F, {{100, -1, 1}}); registry.add("trackTpcNClsCrossedRows", "trackTpcNClsCrossedRows", HistType::kTH1F, {{50, 0, 200}}); @@ -185,6 +206,10 @@ struct JetShapeTask { registry.add("area", "area", HistType::kTH2F, {{100, 0, 2}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("rho", "rho", HistType::kTH2F, {{120, 0, 300}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("ptCorr", "Corrected jet pT; p_{T}^{corr} (GeV/c); Counts", HistType::kTH1F, {{200, 0, 200}}); + registry.add("jetPtCorrCent", "Selected jets;p_{T,jet}^{corr} (GeV/c);Centrality (%)", HistType::kTH2F, {{200, 0, 200}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("areaExposureJet", "Sum of accepted annulus area for signal;r;p_{T,jet}^{corr};Centrality", HistType::kTHnSparseD, {{nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("areaExposurePerp", "Sum of accepted annulus area for two perpendicular cones;r;p_{T,jet}^{corr};Centrality", HistType::kTHnSparseD, {{nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("areaExposureRandCone", "Sum of accepted annulus area for random cones;r;p_{T,jet}^{corr};Centrality", HistType::kTHnSparseD, {{nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("ptCorrVsDistance", "ptcorr_vs_distance", HistType::kTH2F, {{70, 0, 0.7}, {100, 0, 100}}); registry.add("jetDistanceVsTrackpt", "trackpt_vs_distance_injet", HistType::kTH2F, {{70, 0, 0.7}, {100, 0, 100}}); registry.add("ptSum", "ptSum", HistType::kTHnSparseD, {{14, 0, 0.7}, {nBinsJetShapeFunc.value, 0, jetShapeFuncMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); @@ -211,6 +236,10 @@ struct JetShapeTask { registry.add("ptGenPiPerp", "ptGenPiPerp", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("ptGenKaPerp", "ptGenKaPerp", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("ptGenPrPerp", "ptGenPrPerp", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + // Random cone + registry.add("ptGenPiRand", "ptGenPiRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("ptGenKaRand", "ptGenKaRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("ptGenPrRand", "ptGenPrRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); // Inclusive (No JetPt axis) registry.add("ptGenPiInc", "ptGenPiInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("ptGenKaInc", "ptGenKaInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); @@ -224,6 +253,8 @@ struct JetShapeTask { registry.add("effNumTofPiInc", "effNumTofPiInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("effNumTofKaInc", "effNumTofKaInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("effNumTofPrInc", "effNumTofPrInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumPidPiInc", "True primary pion reconstructed with the full analysis PID selection", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumPidPrInc", "True primary proton reconstructed with the full analysis PID selection", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); // Jet (In-cone) registry.add("effNumTpcPiJet", "effNumTpcPiJet", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); @@ -232,6 +263,8 @@ struct JetShapeTask { registry.add("effNumTofPiJet", "effNumTofPiJet", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("effNumTofKaJet", "effNumTofKaJet", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("effNumTofPrJet", "effNumTofPrJet", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumPidPiJet", "Full pion PID efficiency numerator in signal region", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumPidPrJet", "Full proton PID efficiency numerator in signal region", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); // Perp Cone registry.add("effNumTpcPiPerp", "effNumTpcPiPerp", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); @@ -240,6 +273,18 @@ struct JetShapeTask { registry.add("effNumTofPiPerp", "effNumTofPiPerp", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("effNumTofKaPerp", "effNumTofKaPerp", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("effNumTofPrPerp", "effNumTofPrPerp", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumPidPiPerp", "Full pion PID efficiency numerator in perpendicular cones", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumPidPrPerp", "Full proton PID efficiency numerator in perpendicular cones", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + + // Random cone efficiency numerators + registry.add("effNumTpcPiRand", "effNumTpcPiRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumTpcKaRand", "effNumTpcKaRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumTpcPrRand", "effNumTpcPrRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumTofPiRand", "effNumTofPiRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumTofKaRand", "effNumTofKaRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumTofPrRand", "effNumTofPrRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumPidPiRand", "Full pion PID efficiency numerator in random cones", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("effNumPidPrRand", "Full proton PID efficiency numerator in random cones", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); // DCA Templates for Primary Fraction // Jet (In-cone) @@ -258,6 +303,14 @@ struct JetShapeTask { registry.add("dcaDecayPrPerp", "dcaDecayPrPerp", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("dcaMatPrPerp", "dcaMatPrPerp", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + // Random cone + registry.add("dcaPrimPiRand", "dcaPrimPiRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("dcaDecayPiRand", "dcaDecayPiRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("dcaMatPiRand", "dcaMatPiRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("dcaPrimPrRand", "dcaPrimPrRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("dcaDecayPrRand", "dcaDecayPrRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + registry.add("dcaMatPrRand", "dcaMatPrRand", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + // Inclusive registry.add("dcaPrimPiInc", "dcaPrimPiInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("dcaDecayPiInc", "dcaDecayPiInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); @@ -265,6 +318,20 @@ struct JetShapeTask { registry.add("dcaPrimPrInc", "dcaPrimPrInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("dcaDecayPrInc", "dcaDecayPrInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); registry.add("dcaMatPrInc", "dcaMatPrInc", HistType::kTHnSparseD, {{nBinsPt.value, 0, ptMax.value}, {nBinsDcaxyForMc.value, dcaxyMin.value, dcaxyMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}}); + + // Transfer from the wide-DCA fit sample to the nominal global-track sample. + // Axes: + // Jet/Perp: reco pT, reco distance, detector-level corrected jet pT, + // centrality, species (0=pi,1=p), origin (0=primary,1=decay,2=material) + // Inclusive: reco pT, centrality, species, origin + registry.add("originCountWideJet", "wide-DCA origin counts around matched detector-level jets", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}, {2, 0, 2}, {3, 0, 3}}); + registry.add("originCountGlobalJet", "nominal-global origin counts around matched detector-level jets", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}, {2, 0, 2}, {3, 0, 3}}); + registry.add("originCountWidePerp", "wide-DCA origin counts in perpendicular cones", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}, {2, 0, 2}, {3, 0, 3}}); + registry.add("originCountGlobalPerp", "nominal-global origin counts in perpendicular cones", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}, {2, 0, 2}, {3, 0, 3}}); + registry.add("originCountWideRand", "wide-DCA origin counts in random cones", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}, {2, 0, 2}, {3, 0, 3}}); + registry.add("originCountGlobalRand", "nominal-global origin counts in random cones", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsDistance.value, 0, distanceMax.value}, {nBinsJetPt.value, jetPtMinForCut.value, jetPtMaxForCut.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}, {2, 0, 2}, {3, 0, 3}}); + registry.add("originCountWideInc", "inclusive wide-DCA origin counts", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}, {2, 0, 2}, {3, 0, 3}}); + registry.add("originCountGlobalInc", "inclusive nominal-global origin counts", HistType::kTHnSparseD, {{nBinsPtForDca.value, 0, ptMax.value}, {nBinsCentrality.value, centralityMinForCut.value, centralityMaxForCut.value}, {2, 0, 2}, {3, 0, 3}}); } template @@ -281,11 +348,11 @@ struct JetShapeTask { return false; } } - static constexpr double LeadingConstituentPtMinValue = 5.0; static constexpr double LeadingConstituentPtMaxValue = 9998.0; bool checkConstituentPt = true; - bool checkConstituentMinPt = - (leadingConstituentPtMin > LeadingConstituentPtMinValue); + // A non-negative value activates the minimum-leading-constituent cut. + // Therefore the default value 5 GeV/c is actually applied. + bool checkConstituentMinPt = (leadingConstituentPtMin >= 0.0f); bool checkConstituentMaxPt = (leadingConstituentPtMax < LeadingConstituentPtMaxValue); if (!checkConstituentMinPt && !checkConstituentMaxPt) { @@ -312,11 +379,200 @@ struct JetShapeTask { return true; } + template + bool passCommonTrackQuality(T const& track) + { + return track.pt() >= trackPtMinForAnalysis && + std::abs(track.eta()) <= etaTrUp && + track.tpcNClsCrossedRows() >= nclcrossTpcMin && + track.itsChi2NCl() <= chi2ItsMax && + track.tpcChi2NCl() <= chi2TpcMax && + track.tpcNClsFound() >= nclTpcMin && + track.itsNCls() >= nclItsMin; + } + + template + bool passNominalTrackDca(T const& track) + { + return std::abs(track.dcaXY()) <= dcaxyCutMax; + } + + bool passCentrality(float centrality) const + { + return !applyCentralityCut.value || + (centrality >= centralityMinForCut.value && + centrality < centralityMaxForCut.value); + } + + static float absDeltaPhi(float phi1, float phi2) + { + float dPhi = std::abs(phi1 - phi2); + if (dPhi > o2::constants::math::PI) { + dPhi = o2::constants::math::TwoPI - dPhi; + } + return dPhi; + } + + static float rapidityFromPtEtaMass(float pt, float eta, float mass) + { + const float pz = pt * std::sinh(eta); + const float mt = std::sqrt(pt * pt + mass * mass); + return std::asinh(pz / mt); + } + + bool passSpeciesKinematics(float pt, float eta, float mass) const + { + if (std::abs(eta) >= etaTrUp.value) { + return false; + } + + if (!useRapiditySelection.value) { + return true; + } + + return std::abs(rapidityFromPtEtaMass(pt, eta, mass)) < + mcRapidityMax.value; + } + + template + bool passTruthKinematics(T const& particle) const + { + if (particle.pt() < trackPtMinForAnalysis.value || + std::abs(particle.eta()) >= etaTrUp.value) { + return false; + } + + return !useRapiditySelection.value || + std::abs(particle.y()) < mcRapidityMax.value; + } + + static double diskAreaPrimitive(double radius, double x) + { + if (radius <= 0.0) { + return 0.0; + } + x = std::clamp(x, -radius, radius); + const double root = std::sqrt(std::max(0.0, radius * radius - x * x)); + return x * root + radius * radius * std::asin(x / radius); + } + + static double effectiveDiskArea(double radius, double axisEta, double etaMin, double etaMax) + { + const double lower = std::max(-radius, etaMin - axisEta); + const double upper = std::min(radius, etaMax - axisEta); + if (upper <= lower) { + return 0.0; + } + return diskAreaPrimitive(radius, upper) - diskAreaPrimitive(radius, lower); + } + + static double effectiveAnnulusArea(double rMin, double rMax, double axisEta, double etaMin, double etaMax) + { + return effectiveDiskArea(rMax, axisEta, etaMin, etaMax) - + effectiveDiskArea(rMin, axisEta, etaMin, etaMax); + } + + struct ConeVetoJet { + float eta; + float phi; + float ptCorr; + }; + + struct RandomConeAxis { + float eta; + float phi; + float referenceJetPtCorr; + }; + + std::vector makeRandomCones(std::vector const& jets) + { + std::vector cones; + if (!enableRandomCone || jets.empty() || etaTrUp <= distanceMax) { + return cones; + } + + cones.reserve(jets.size()); + const float etaMin = -etaTrUp + distanceMax; + const float etaMax = etaTrUp - distanceMax; + const float vetoDistance = jetR + distanceMax; + + for (const auto& referenceJet : jets) { + bool accepted = false; + float rcEta = 0.0f; + float rcPhi = 0.0f; + + for (int attempt = 0; attempt < randomConeMaxAttempts; ++attempt) { + rcEta = randomNumber.Uniform(etaMin, etaMax); + if (randomConeFullAzimuth) { + rcPhi = randomNumber.Uniform(0.0f, static_cast(o2::constants::math::TwoPI)); + } else { + float dPhi = randomNumber.Uniform(randomConeDeltaPhiMin, randomConeDeltaPhiMax); + if (randomNumber.Uniform() < RandomConeNegativeSideProbability) { + dPhi = -dPhi; + } + rcPhi = RecoDecay::constrainAngle(referenceJet.phi + dPhi); + } + + accepted = true; + for (const auto& vetoJet : jets) { + const float dEta = rcEta - vetoJet.eta; + const float dPhi = absDeltaPhi(rcPhi, vetoJet.phi); + if (dEta * dEta + dPhi * dPhi <= vetoDistance * vetoDistance) { + accepted = false; + break; + } + } + if (accepted) { + break; + } + } + + if (accepted) { + cones.push_back({rcEta, rcPhi, referenceJet.ptCorr}); + } + } + return cones; + } + + enum class ExposureRegion { Signal, + Perpendicular }; + + void fillAreaExposure(float axisEta, float jetPtCorr, float centrality, ExposureRegion region, double multiplicity = 1.0) + { + const double dr = distanceMax / static_cast(nBinsDistance); + for (int i = 0; i < nBinsDistance; ++i) { + const double rMin = i * dr; + const double rMax = (i + 1) * dr; + const double rCenter = 0.5 * (rMin + rMax); + const double area = multiplicity * effectiveAnnulusArea(rMin, rMax, axisEta, -etaTrUp, etaTrUp); + if (region == ExposureRegion::Signal) { + registry.fill(HIST("areaExposureJet"), rCenter, jetPtCorr, centrality, area); + } else { + registry.fill(HIST("areaExposurePerp"), rCenter, jetPtCorr, centrality, area); + } + } + } + + void fillFullRandomConeArea(float jetPtCorr, float centrality) + { + const double dr = distanceMax / static_cast(nBinsDistance); + for (int i = 0; i < nBinsDistance; ++i) { + const double rMin = i * dr; + const double rMax = (i + 1) * dr; + const double rCenter = 0.5 * (rMin + rMax); + const double area = o2::constants::math::PI * (rMax * rMax - rMin * rMin); + registry.fill(HIST("areaExposureRandCone"), rCenter, jetPtCorr, centrality, area); + } + } + Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); Filter jetCollisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; Filter collisionFilter = nabs(aod::collision::posZ) < vertexZCut; Filter mcCollisionFilter = nabs(aod::jmccollision::posZ) < vertexZCut; + using JetTracksWithPI = soa::Join; + using JetTracksWithPIMC = soa::Join; + using FullTrackInfo = soa::Join; using FullTrackInfoMC = soa::Join; @@ -326,7 +582,11 @@ struct JetShapeTask { return; } - registry.fill(HIST("eventCounterJetShape"), collision.centFT0M()); + const float cent = collision.centFT0M(); + if (!passCentrality(cent)) { + return; + } + registry.fill(HIST("eventCounterJetShape"), cent); size_t nBins = distanceCategory->size() - 1; @@ -349,7 +609,6 @@ struct JetShapeTask { accumSumBg2.reserve(jets.size()); float rho = collision.rho(); - float cent = collision.centFT0M(); for (const auto& jet : jets) { if (!isAcceptedJet(jet)) { @@ -371,6 +630,10 @@ struct JetShapeTask { accumSumBg2.push_back(std::vector(nBins, 0.f)); } + if (cachedJets.empty()) { + return; + } + for (const auto& track : tracks) { if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { continue; @@ -385,8 +648,9 @@ struct JetShapeTask { float dEta = trkEta - jet.eta; - if (std::abs(dEta) > maxDistance) + if (std::abs(dEta) > maxDistance) { continue; + } float dPhi = RecoDecay::constrainAngle(trkPhi - jet.phi); float distance = std::sqrt(dEta * dEta + dPhi * dPhi); @@ -465,318 +729,341 @@ struct JetShapeTask { PROCESS_SWITCH(JetShapeTask, processJetShape, "JetShape", false); - void processJetProductionRatio(soa::Filtered>::iterator const& collision, soa::Join const& tracks, FullTrackInfo const&, soa::Join const& jets) + template + void runJetProductionRatio(CollisionT const& collision, + JetTracksWithPI const& tracks, + soa::Join const& jets, + float rho) { if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; } - registry.fill(HIST("eventCounterJet"), collision.centFT0M()); - registry.fill(HIST("event/vertexz"), collision.posZ()); - float rho = collision.rho(); - float centrality = collision.centFT0M(); + const float centrality = collision.centFT0M(); + if (!passCentrality(centrality)) { + return; + } + + registry.fill(HIST("eventCounterJet"), centrality); + registry.fill(HIST("event/vertexz"), collision.posZ()); + registry.fill(HIST("rho"), rho, centrality); struct CachedJet { float pt, eta, phi, ptCorr; float phiBg1, phiBg2; }; std::vector cachedJets; + std::vector vetoJets; cachedJets.reserve(jets.size()); + vetoJets.reserve(jets.size()); for (const auto& jet : jets) { - if (!isAcceptedJet(jet)) + if (!isAcceptedJet(jet)) { continue; + } - float ptCorr = jet.pt() - rho * jet.area(); - float phiBg1 = - RecoDecay::constrainAngle(jet.phi() + o2::constants::math::PIHalf); - float phiBg2 = - RecoDecay::constrainAngle(jet.phi() - o2::constants::math::PIHalf); + const float ptCorr = jet.pt() - rho * jet.area(); + const float phiBg1 = RecoDecay::constrainAngle(jet.phi() + o2::constants::math::PIHalf); + const float phiBg2 = RecoDecay::constrainAngle(jet.phi() - o2::constants::math::PIHalf); registry.fill(HIST("jetEta"), jet.eta()); registry.fill(HIST("jetPhi"), jet.phi()); registry.fill(HIST("area"), jet.area(), centrality); - registry.fill(HIST("rho"), rho, centrality); registry.fill(HIST("jetPt"), jet.pt(), centrality); + registry.fill(HIST("ptCorr"), ptCorr); + registry.fill(HIST("jetPtCorrCent"), ptCorr, centrality); - cachedJets.push_back( - {jet.pt(), jet.eta(), jet.phi(), ptCorr, phiBg1, phiBg2}); - } - - bool isValidRC = false; - float rcEta = 0.0f; - float rcPhi = 0.0f; - - if (!cachedJets.empty()) { - constexpr unsigned int RandomSeed = 0; - TRandom3 randomNumber(RandomSeed); - - const auto& leadJet = cachedJets[0]; - - constexpr float MaxTrackEta = 0.9f; - constexpr int MaxAttempts = 100; // for RandomCone - constexpr float FlipProbability = 0.5f; - - // Range to generate - float rcEtaMin = -MaxTrackEta + distanceMax; - float rcEtaMax = MaxTrackEta - distanceMax; - - int attempts = 0; - while (!isValidRC && attempts < MaxAttempts) { - rcEta = randomNumber.Uniform(rcEtaMin, rcEtaMax); - - float dPhi = randomNumber.Uniform(randomConeDeltaPhiMin, randomConeDeltaPhiMax); - - // flipProbability (0.5) - if (randomNumber.Uniform() < FlipProbability) { - dPhi = -dPhi; - } - - rcPhi = RecoDecay::constrainAngle(leadJet.phi + dPhi); + fillAreaExposure(jet.eta(), ptCorr, centrality, ExposureRegion::Signal); + fillAreaExposure(jet.eta(), ptCorr, centrality, ExposureRegion::Perpendicular, 2.0); - float dPhiLead = std::abs(rcPhi - leadJet.phi); - if (dPhiLead > o2::constants::math::PI) - dPhiLead = o2::constants::math::TwoPI - dPhiLead; - float dEtaLead = rcEta - leadJet.eta; - float distLead = std::sqrt(dEtaLead * dEtaLead + dPhiLead * dPhiLead); + cachedJets.push_back({jet.pt(), jet.eta(), jet.phi(), ptCorr, phiBg1, phiBg2}); + vetoJets.push_back({jet.eta(), jet.phi(), ptCorr}); + } - if (distLead > (jetR + distanceMax)) { - isValidRC = true; - } - attempts++; - } + if (cachedJets.empty()) { + return; + } - if (isValidRC) { - const auto& leadJet = cachedJets[0]; - registry.fill(HIST("eventCounterRandCone"), leadJet.ptCorr, centrality); - } + const auto randomCones = makeRandomCones(vetoJets); + for (const auto& cone : randomCones) { + registry.fill(HIST("eventCounterRandCone"), cone.referenceJetPtCorr, centrality); + fillFullRandomConeArea(cone.referenceJetPtCorr, centrality); } + const float maxR2 = distanceMax * distanceMax; + + // The jet finder remains configured with globalTracks. The nominal-yield + // sample and the wide-DCA sample are selected independently relative to + // the same reconstructed jet axes. for (const auto& jetTrack : tracks) { - if (!jetderiveddatautilities::selectTrack(jetTrack, trackSelection)) { + const bool passesNominalBit = jetderiveddatautilities::selectTrack(jetTrack, trackSelection); + const bool passesWideDcaBit = jetderiveddatautilities::selectTrack(jetTrack, trackSelectionForDca); + if (!passesNominalBit && !passesWideDcaBit) { continue; } - auto track = jetTrack.track_as(); - - if (std::abs(track.eta()) > etaTrUp) - continue; - if (track.tpcNClsCrossedRows() < nclcrossTpcMin) - continue; - if (std::abs(track.dcaXY()) > dcaxyCutMax) - continue; - if (track.itsChi2NCl() > chi2ItsMax) + const auto track = jetTrack.track_as(); + if (!passCommonTrackQuality(track)) { continue; - if (track.tpcChi2NCl() > chi2TpcMax) - continue; - if (track.tpcNClsFound() < nclTpcMin) - continue; - if (track.itsNCls() < nclItsMin) - continue; - - float trkPt = track.pt(); - float trkP = track.p(); - float trkEta = track.eta(); - float trkPhi = track.phi(); - - float tpcSig = track.tpcSignal(); - float tpcPi = track.tpcNSigmaPi(); - float tofPi = track.tofNSigmaPi(); - float tpcPr = track.tpcNSigmaPr(); - float tofPr = track.tofNSigmaPr(); - float beta = track.beta(); + } - // PID QA - registry.fill(HIST("jetTpcPi"), trkP, tpcPi); - registry.fill(HIST("jetTofPi"), trkPt, tofPi); - registry.fill(HIST("jetTpcPr"), trkP, tpcPr); - registry.fill(HIST("jetTofPr"), trkPt, tofPr); + const bool passesNominal = passesNominalBit && passNominalTrackDca(track); + const bool passesWideDca = passesWideDcaBit; + if (!passesNominal && !passesWideDca) { + continue; + } - bool hasTofPi = (std::abs(tofPi) < nSigmaTofCut); - bool hasTofPr = (std::abs(tofPr) < nSigmaTofCut); - bool isTpcPiRange = (tpcPi > tpcNSigmaPiMin && tpcPi < tpcNSigmaPiMax); - bool isTpcPrRange = (tpcPr > tpcNSigmaPrMin && tpcPr < tpcNSigmaPrMax); + const float trkPt = track.pt(); + const float trkP = track.p(); + const float trkEta = track.eta(); + const float trkPhi = track.phi(); + + const float tpcSig = track.tpcSignal(); + const float tpcPi = track.tpcNSigmaPi(); + const float tofPi = track.tofNSigmaPi(); + const float tpcPr = track.tpcNSigmaPr(); + const float tofPr = track.tofNSigmaPr(); + const float beta = track.beta(); + + const bool passPiKinematics = passSpeciesKinematics(trkPt, trkEta, PionMass); + const bool passPrKinematics = passSpeciesKinematics(trkPt, trkEta, ProtonMass); + const bool hasTofPi = passPiKinematics && track.hasTOF() && (std::abs(tofPi) < nSigmaTofCut); + const bool hasTofPr = passPrKinematics && track.hasTOF() && (std::abs(tofPr) < nSigmaTofCut); + const bool isTpcPiRange = (tpcPi > tpcNSigmaPiMin && tpcPi < tpcNSigmaPiMax); + const bool isTpcPrRange = (tpcPr > tpcNSigmaPrMin && tpcPr < tpcNSigmaPrMax); + + if (passesNominal && fillDetailedQA) { + registry.fill(HIST("jetTpcPi"), trkP, tpcPi); + registry.fill(HIST("jetTofPi"), trkPt, tofPi); + registry.fill(HIST("jetTpcPr"), trkP, tpcPr); + registry.fill(HIST("jetTofPr"), trkPt, tofPr); + } + // Signal and perpendicular cones: one loop over accepted jets. for (const auto& jet : cachedJets) { - - float dEta = trkEta - jet.eta; - if (std::abs(dEta) > maxDeltaEtaSafe) + const float dEta = trkEta - jet.eta; + if (std::abs(dEta) >= distanceMax) { continue; - - float dPhi = RecoDecay::constrainAngle(trkPhi - jet.phi); - float distance = std::sqrt(dEta * dEta + dPhi * dPhi); - - // Background judge - float deltaPhiBg1 = RecoDecay::constrainAngle(trkPhi - jet.phiBg1); - float deltaPhiBg2 = RecoDecay::constrainAngle(trkPhi - jet.phiBg2); - float distBg1 = std::sqrt(dEta * dEta + deltaPhiBg1 * deltaPhiBg1); - float distBg2 = std::sqrt(dEta * dEta + deltaPhiBg2 * deltaPhiBg2); - - float distBg = -1.0f; - if (distBg1 < distanceMax) { - distBg = distBg1; - } else if (distBg2 < distanceMax) { - distBg = distBg2; } - // --- Background Fill --- - if (distBg >= 0.0f) { - registry.fill(HIST("tpcDedxPerpJet"), trkP, tpcSig); + const float dPhi = absDeltaPhi(trkPhi, jet.phi); + const float distanceSq = dEta * dEta + dPhi * dPhi; - // dcaXY - if (track.hasTOF()) { + auto fillPerpendicular = [&](float conePhi) { + const float dPhiBg = absDeltaPhi(trkPhi, conePhi); + const float distBgSq = dEta * dEta + dPhiBg * dPhiBg; + if (distBgSq >= maxR2) { + return; + } + const float distBg = std::sqrt(distBgSq); + + if (passesWideDca) { if (hasTofPr && isTpcPrRange) { registry.fill(HIST("dcaPrPerpJet"), trkPt, track.dcaXY(), distBg, jet.ptCorr, centrality); } - if (hasTofPi && isTpcPiRange) { registry.fill(HIST("dcaPiPerpJet"), trkPt, track.dcaXY(), distBg, jet.ptCorr, centrality); } } - if (hasTofPi) { - registry.fill(HIST("tpcTofPiPerpJet"), trkP, tpcPi, distBg, jet.ptCorr, centrality); - if (isTpcPiRange) { - registry.fill(HIST("pVsPtForPiPerpJet"), trkP, trkPt, distBg, jet.ptCorr, centrality); - } - } - if (hasTofPr) { - registry.fill(HIST("tpcTofPrPerpJet"), trkP, tpcPr, distBg, jet.ptCorr, centrality); - if (isTpcPrRange) { - registry.fill(HIST("pVsPtForPrPerpJet"), trkP, trkPt, distBg, jet.ptCorr, centrality); - } - } - } - - if (isValidRC) { - const auto& leadJet = cachedJets[0]; - - float dEtaRC = trkEta - rcEta; - float dPhiRC = std::abs(trkPhi - rcPhi); - if (dPhiRC > o2::constants::math::PI) - dPhiRC = o2::constants::math::TwoPI - dPhiRC; - float distRC = std::sqrt(dEtaRC * dEtaRC + dPhiRC * dPhiRC); - - if (distRC < distanceMax) { - - // dcaXY - if (track.hasTOF()) { - if (hasTofPr && isTpcPrRange) { - registry.fill(HIST("dcaPrRandCone"), trkPt, track.dcaXY(), distRC, leadJet.ptCorr, centrality); - } - if (hasTofPi && isTpcPiRange) { - registry.fill(HIST("dcaPiRandCone"), trkPt, track.dcaXY(), distRC, leadJet.ptCorr, centrality); - } + if (passesNominal) { + if (fillDetailedQA) { + registry.fill(HIST("tpcDedxPerpJet"), trkP, tpcSig); } - if (hasTofPi) { - registry.fill(HIST("tpcTofPiRandCone"), trkP, tpcPi, distRC, leadJet.ptCorr, centrality); + registry.fill(HIST("tpcTofPiPerpJet"), trkP, tpcPi, distBg, jet.ptCorr, centrality); if (isTpcPiRange) { - registry.fill(HIST("pVsPtForPiRandCone"), trkP, trkPt, distRC, leadJet.ptCorr, centrality); + registry.fill(HIST("pVsPtForPiPerpJet"), trkP, trkPt, distBg, jet.ptCorr, centrality); } } if (hasTofPr) { - registry.fill(HIST("tpcTofPrRandCone"), trkP, tpcPr, distRC, leadJet.ptCorr, centrality); + registry.fill(HIST("tpcTofPrPerpJet"), trkP, tpcPr, distBg, jet.ptCorr, centrality); if (isTpcPrRange) { - registry.fill(HIST("pVsPtForPrRandCone"), trkP, trkPt, distRC, leadJet.ptCorr, centrality); + registry.fill(HIST("pVsPtForPrPerpJet"), trkP, trkPt, distBg, jet.ptCorr, centrality); } } } - } + }; + + fillPerpendicular(jet.phiBg1); + fillPerpendicular(jet.phiBg2); - // --- Signal Fill --- - registry.fill(HIST("jetDistanceVsTrackpt"), distance, trkPt); - registry.fill(HIST("jetTpcDedx"), trkP, tpcSig, distance); - registry.fill(HIST("jetTofBeta"), trkP, beta); + if (distanceSq >= maxR2) { + continue; + } + const float distance = std::sqrt(distanceSq); - // dcaXY - if (track.hasTOF()) { + if (passesWideDca) { if (hasTofPr && isTpcPrRange) { registry.fill(HIST("jetDcaPr"), trkPt, track.dcaXY(), distance, jet.ptCorr, centrality); } - if (hasTofPi && isTpcPiRange) { registry.fill(HIST("jetDcaPi"), trkPt, track.dcaXY(), distance, jet.ptCorr, centrality); } } - if (hasTofPr) { - registry.fill(HIST("jetTpcTofPr"), trkP, tpcPr, distance, jet.ptCorr, centrality); - if (isTpcPrRange) { - registry.fill(HIST("jetpVsPtForPr"), trkP, trkPt, distance, jet.ptCorr, centrality); + if (passesNominal) { + if (fillDetailedQA) { + registry.fill(HIST("jetDistanceVsTrackpt"), distance, trkPt); + registry.fill(HIST("jetTpcDedx"), trkP, tpcSig, distance); + registry.fill(HIST("jetTofBeta"), trkP, beta); + } + + if (hasTofPr) { + registry.fill(HIST("jetTpcTofPr"), trkP, tpcPr, distance, jet.ptCorr, centrality); + if (isTpcPrRange) { + registry.fill(HIST("jetpVsPtForPr"), trkP, trkPt, distance, jet.ptCorr, centrality); + } + } + if (hasTofPi) { + registry.fill(HIST("jetTpcTofPi"), trkP, tpcPi, distance, jet.ptCorr, centrality); + if (isTpcPiRange) { + registry.fill(HIST("jetpVsPtForPi"), trkP, trkPt, distance, jet.ptCorr, centrality); + } + } + } + } + + // Random cones are outside the jet loop, so every track-cone pair is + // filled exactly once. Each cone is associated with one selected jet. + for (const auto& cone : randomCones) { + const float dEtaRC = trkEta - cone.eta; + if (std::abs(dEtaRC) >= distanceMax) { + continue; + } + const float dPhiRC = absDeltaPhi(trkPhi, cone.phi); + const float distRCSq = dEtaRC * dEtaRC + dPhiRC * dPhiRC; + if (distRCSq >= maxR2) { + continue; + } + const float distRC = std::sqrt(distRCSq); + + if (passesWideDca) { + if (hasTofPr && isTpcPrRange) { + registry.fill(HIST("dcaPrRandCone"), trkPt, track.dcaXY(), distRC, cone.referenceJetPtCorr, centrality); + } + if (hasTofPi && isTpcPiRange) { + registry.fill(HIST("dcaPiRandCone"), trkPt, track.dcaXY(), distRC, cone.referenceJetPtCorr, centrality); } } - if (hasTofPi) { - registry.fill(HIST("jetTpcTofPi"), trkP, tpcPi, distance, jet.ptCorr, centrality); - if (isTpcPiRange) { - registry.fill(HIST("jetpVsPtForPi"), trkP, trkPt, distance, jet.ptCorr, centrality); + if (passesNominal) { + if (hasTofPi) { + registry.fill(HIST("tpcTofPiRandCone"), trkP, tpcPi, distRC, cone.referenceJetPtCorr, centrality); + if (isTpcPiRange) { + registry.fill(HIST("pVsPtForPiRandCone"), trkP, trkPt, distRC, cone.referenceJetPtCorr, centrality); + } + } + if (hasTofPr) { + registry.fill(HIST("tpcTofPrRandCone"), trkP, tpcPr, distRC, cone.referenceJetPtCorr, centrality); + if (isTpcPrRange) { + registry.fill(HIST("pVsPtForPrRandCone"), trkP, trkPt, distRC, cone.referenceJetPtCorr, centrality); + } } } } } } - PROCESS_SWITCH(JetShapeTask, processJetProductionRatio, "production ratio around jets", false); - void processInclusiveProductionRatio(soa::Filtered::iterator const& collision, soa::Join const& tracks, FullTrackInfo const&) + void processJetProductionRatioWithRho( + soa::Filtered>::iterator const& collision, + JetTracksWithPI const& tracks, + FullTrackInfo const&, + soa::Join const& jets) + { + runJetProductionRatio(collision, tracks, jets, collision.rho()); + } + PROCESS_SWITCH(JetShapeTask, processJetProductionRatioWithRho, "production ratio around jets with rho subtraction", false); + + void processJetProductionRatioNoRho( + soa::Filtered::iterator const& collision, + JetTracksWithPI const& tracks, + FullTrackInfo const&, + soa::Join const& jets) + { + runJetProductionRatio(collision, tracks, jets, 0.0f); + } + PROCESS_SWITCH(JetShapeTask, processJetProductionRatioNoRho, "production ratio around jets without rho table (recommended for pp)", false); + + void processInclusiveProductionRatio(soa::Filtered::iterator const& collision, JetTracksWithPI const& tracks, FullTrackInfo const&) { if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; } - registry.fill(HIST("eventCounterInc"), collision.centFT0M()); - // tracks conditions + const float centrality = collision.centFT0M(); + if (!passCentrality(centrality)) { + return; + } + registry.fill(HIST("eventCounterInc"), centrality); + for (const auto& jetTrack : tracks) { + const bool passesNominalBit = jetderiveddatautilities::selectTrack(jetTrack, trackSelection); + const bool passesWideDcaBit = jetderiveddatautilities::selectTrack(jetTrack, trackSelectionForDca); + if (!passesNominalBit && !passesWideDcaBit) { + continue; + } - if (!jetderiveddatautilities::selectTrack(jetTrack, trackSelection)) { + const auto track = jetTrack.track_as(); + + // Keep the original QA tied to the nominal global-track selection. + if (passesNominalBit && fillDetailedQA) { + registry.fill(HIST("trackTpcNClsCrossedRows"), track.tpcNClsCrossedRows()); + registry.fill(HIST("trackDcaXY"), track.dcaXY()); + registry.fill(HIST("trackItsChi2NCl"), track.itsChi2NCl()); + registry.fill(HIST("trackTpcChi2NCl"), track.tpcChi2NCl()); + registry.fill(HIST("trackTpcNClsFound"), track.tpcNClsFound()); + registry.fill(HIST("trackItsNCls"), track.itsNCls()); + registry.fill(HIST("trackEta"), track.eta()); + registry.fill(HIST("trackPhi"), track.phi()); + } + + if (!passCommonTrackQuality(track)) { continue; } - auto track = jetTrack.track_as(); + const bool passesNominal = passesNominalBit && passNominalTrackDca(track); + const bool passesWideDca = passesWideDcaBit; - registry.fill(HIST("trackTpcNClsCrossedRows"), track.tpcNClsCrossedRows()); - registry.fill(HIST("trackDcaXY"), track.dcaXY()); - registry.fill(HIST("trackItsChi2NCl"), track.itsChi2NCl()); - registry.fill(HIST("trackTpcChi2NCl"), track.tpcChi2NCl()); - registry.fill(HIST("trackTpcNClsFound"), track.tpcNClsFound()); - registry.fill(HIST("trackItsNCls"), track.itsNCls()); - registry.fill(HIST("trackEta"), track.eta()); - registry.fill(HIST("trackPhi"), track.phi()); + const bool passPiKinematics = passSpeciesKinematics(track.pt(), track.eta(), PionMass); + const bool passPrKinematics = passSpeciesKinematics(track.pt(), track.eta(), ProtonMass); + const bool hasTofPi = passPiKinematics && track.hasTOF() && std::abs(track.tofNSigmaPi()) < nSigmaTofCut; + const bool hasTofPr = passPrKinematics && track.hasTOF() && std::abs(track.tofNSigmaPr()) < nSigmaTofCut; + const bool isTpcPiRange = track.tpcNSigmaPi() > tpcNSigmaPiMin && track.tpcNSigmaPi() < tpcNSigmaPiMax; + const bool isTpcPrRange = track.tpcNSigmaPr() > tpcNSigmaPrMin && track.tpcNSigmaPr() < tpcNSigmaPrMax; - if (std::abs(track.eta()) > etaTrUp) - continue; - if (track.tpcNClsCrossedRows() < nclcrossTpcMin) - continue; - if (std::abs(track.dcaXY()) > dcaxyCutMax) - continue; - if (track.itsChi2NCl() > chi2ItsMax) - continue; - if (track.tpcChi2NCl() > chi2TpcMax) - continue; - if (track.tpcNClsFound() < nclTpcMin) - continue; - if (track.itsNCls() < nclItsMin) - continue; + if (passesWideDca) { + if (hasTofPi && isTpcPiRange) { + registry.fill(HIST("dcaPiInc"), track.pt(), track.dcaXY(), centrality); + } + if (hasTofPr && isTpcPrRange) { + registry.fill(HIST("dcaPrInc"), track.pt(), track.dcaXY(), centrality); + } + } - // PID check - registry.fill(HIST("tpcPi"), track.p(), track.tpcNSigmaPi()); - registry.fill(HIST("tofPi"), track.pt(), track.tofNSigmaPi()); - registry.fill(HIST("tpcPr"), track.p(), track.tpcNSigmaPr()); - registry.fill(HIST("tofPr"), track.pt(), track.tofNSigmaPr()); - registry.fill(HIST("tpcDedx"), track.p(), track.tpcSignal(), collision.centFT0M()); - registry.fill(HIST("tofBeta"), track.p(), track.beta(), collision.centFT0M()); + if (!passesNominal) { + continue; + } - if (std::abs(track.tofNSigmaPr()) < nSigmaTofCut) { - registry.fill(HIST("tpcTofPr"), track.p(), track.tpcNSigmaPr(), collision.centFT0M()); + if (fillDetailedQA) { + registry.fill(HIST("tpcPi"), track.p(), track.tpcNSigmaPi()); + registry.fill(HIST("tofPi"), track.pt(), track.tofNSigmaPi()); + registry.fill(HIST("tpcPr"), track.p(), track.tpcNSigmaPr()); + registry.fill(HIST("tofPr"), track.pt(), track.tofNSigmaPr()); + registry.fill(HIST("tpcDedx"), track.p(), track.tpcSignal(), centrality); + registry.fill(HIST("tofBeta"), track.p(), track.beta(), centrality); + } - if (track.tpcNSigmaPr() > tpcNSigmaPrMin && track.tpcNSigmaPr() < tpcNSigmaPrMax) { - registry.fill(HIST("pVsPtForPr"), track.p(), track.pt(), collision.centFT0M()); + if (hasTofPr) { + registry.fill(HIST("tpcTofPr"), track.p(), track.tpcNSigmaPr(), centrality); + if (isTpcPrRange) { + registry.fill(HIST("pVsPtForPr"), track.p(), track.pt(), centrality); } } - if (std::abs(track.tofNSigmaPi()) < nSigmaTofCut) { - registry.fill(HIST("tpcTofPi"), track.p(), track.tpcNSigmaPi(), collision.centFT0M()); - if (track.tpcNSigmaPi() > tpcNSigmaPiMin && track.tpcNSigmaPi() < tpcNSigmaPiMax) { - registry.fill(HIST("pVsPtForPi"), track.p(), track.pt(), collision.centFT0M()); + if (hasTofPi) { + registry.fill(HIST("tpcTofPi"), track.p(), track.tpcNSigmaPi(), centrality); + if (isTpcPiRange) { + registry.fill(HIST("pVsPtForPi"), track.p(), track.pt(), centrality); } } } @@ -785,321 +1072,542 @@ struct JetShapeTask { void processEfficiencyAndPurity( soa::Filtered>::iterator const& collision, - soa::Join const& jetTracks, + JetTracksWithPIMC const& jetTracks, FullTrackInfoMC const&, soa::Filtered> const& mcdJets, aod::ChargedMCParticleLevelJets const& mcpJets, - aod::McParticles const& mcParticles) + aod::JetParticles const& mcParticles) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + return; + } + if (!collision.has_mcCollision()) { return; + } + const auto reducedMcCollisionId = collision.mcCollisionId(); (void)mcpJets; - float centrality = collision.centFT0M(); - float rho = collision.rho(); - float distMax = distanceMax; - float maxR2 = distMax * distMax; + const float centrality = collision.centFT0M(); + if (!passCentrality(centrality)) { + return; + } + const float rho = collision.rho(); + const float maxR2 = distanceMax * distanceMax; registry.fill(HIST("eventCounterMc"), centrality); - // jet matching struct MatchedJet { - float detPtCorr, partEta, partPhi, partPhiBg1, partPhiBg2; + float detPtCorr; + float detEta; + float detPhi; + float detPhiBg1; + float detPhiBg2; + float partEta; + float partPhi; + float partPhiBg1; + float partPhiBg2; }; + std::vector validJets; + std::vector vetoJets; validJets.reserve(mcdJets.size()); + vetoJets.reserve(mcdJets.size()); - for (auto const& detJet : mcdJets) { - auto matchedIndices = detJet.matchedJetGeo(); - - if (matchedIndices.size() > 0) { - auto const& partJet = matchedIndices[0]; + for (const auto& detJet : mcdJets) { + if (!isAcceptedJet(detJet)) { + continue; + } - float detPtCorr = detJet.pt() - rho * detJet.area(); - float partPhi = partJet.phi(); - float phiBg1 = RecoDecay::constrainAngle(partPhi + o2::constants::math::PIHalf); - float phiBg2 = RecoDecay::constrainAngle(partPhi - o2::constants::math::PIHalf); + const auto matchedJets = detJet.matchedJetGeo(); + if (matchedJets.empty()) { + continue; + } - validJets.push_back({detPtCorr, partJet.eta(), partPhi, phiBg1, phiBg2}); - registry.fill(HIST("jetPtMc"), detJet.pt(), centrality); + const auto& partJet = matchedJets[0]; + if (partJet.mcCollisionId() != reducedMcCollisionId) { + continue; } + + const float detPtCorr = detJet.pt() - rho * detJet.area(); + validJets.push_back({detPtCorr, + detJet.eta(), + detJet.phi(), + RecoDecay::constrainAngle(detJet.phi() + o2::constants::math::PIHalf), + RecoDecay::constrainAngle(detJet.phi() - o2::constants::math::PIHalf), + partJet.eta(), + partJet.phi(), + RecoDecay::constrainAngle(partJet.phi() + o2::constants::math::PIHalf), + RecoDecay::constrainAngle(partJet.phi() - o2::constants::math::PIHalf)}); + + vetoJets.push_back({detJet.eta(), detJet.phi(), detPtCorr}); + registry.fill(HIST("jetPtMc"), detJet.pt(), centrality); } - // Denominator: True Primary Particles - for (const auto& mcParticle : mcParticles) { + const auto randomCones = makeRandomCones(vetoJets); - if (mcParticle.mcCollisionId() != collision.globalIndex()) + // Generated denominator for tracking efficiency + for (const auto& mcParticle : mcParticles) { + if (mcParticle.mcCollisionId() != reducedMcCollisionId) { continue; - - if (!mcParticle.isPhysicalPrimary() || std::abs(mcParticle.y()) > mcRapidityMax) + } + if (!mcParticle.isPhysicalPrimary()) { continue; + } + if (!passTruthKinematics(mcParticle)) { + continue; + } - int absPdg = std::abs(mcParticle.pdgCode()); - bool isPi = (absPdg == PDG_t::kPiPlus); - bool isKa = (absPdg == PDG_t::kKPlus); - bool isPr = (absPdg == PDG_t::kProton); - if (!isPi && !isKa && !isPr) + const int absPdg = std::abs(mcParticle.pdgCode()); + const bool isPi = absPdg == PDG_t::kPiPlus; + const bool isKa = absPdg == PDG_t::kKPlus; + const bool isPr = absPdg == PDG_t::kProton; + if (!isPi && !isKa && !isPr) { continue; + } - float mcPt = mcParticle.pt(); + const float mcPt = mcParticle.pt(); - if (isPi) + if (isPi) { registry.fill(HIST("ptGenPiInc"), mcPt, centrality); - else if (isKa) + } else if (isKa) { registry.fill(HIST("ptGenKaInc"), mcPt, centrality); - else if (isPr) + } else { registry.fill(HIST("ptGenPrInc"), mcPt, centrality); + } for (const auto& jet : validJets) { - float dEta = mcParticle.eta() - jet.partEta; - - // Jet In-cone - float dPhiJet = std::abs(mcParticle.phi() - jet.partPhi); - if (dPhiJet > o2::constants::math::PI) - dPhiJet = o2::constants::math::TwoPI - dPhiJet; - - float distJetGen = std::sqrt(dEta * dEta + dPhiJet * dPhiJet); - - if ((dEta * dEta + dPhiJet * dPhiJet) < maxR2) { - if (isPi) - registry.fill(HIST("ptGenPi"), mcPt, distJetGen, jet.detPtCorr, centrality); - else if (isKa) - registry.fill(HIST("ptGenKa"), mcPt, distJetGen, jet.detPtCorr, centrality); - else if (isPr) - registry.fill(HIST("ptGenPr"), mcPt, distJetGen, jet.detPtCorr, centrality); + const float dEtaPart = mcParticle.eta() - jet.partEta; + + const float dPhiJetPart = absDeltaPhi(mcParticle.phi(), jet.partPhi); + const float distJetSq = dEtaPart * dEtaPart + dPhiJetPart * dPhiJetPart; + if (distJetSq < maxR2) { + const float distJet = std::sqrt(distJetSq); + if (isPi) { + registry.fill(HIST("ptGenPi"), mcPt, distJet, jet.detPtCorr, centrality); + } else if (isKa) { + registry.fill(HIST("ptGenKa"), mcPt, distJet, jet.detPtCorr, centrality); + } else { + registry.fill(HIST("ptGenPr"), mcPt, distJet, jet.detPtCorr, centrality); + } + } + + const float dPhiBg1Part = absDeltaPhi(mcParticle.phi(), jet.partPhiBg1); + const float dPhiBg2Part = absDeltaPhi(mcParticle.phi(), jet.partPhiBg2); + const float distBg1Sq = dEtaPart * dEtaPart + dPhiBg1Part * dPhiBg1Part; + const float distBg2Sq = dEtaPart * dEtaPart + dPhiBg2Part * dPhiBg2Part; + + float distBg = -1.0f; + if (distBg1Sq < maxR2) { + distBg = std::sqrt(distBg1Sq); + } else if (distBg2Sq < maxR2) { + distBg = std::sqrt(distBg2Sq); + } + + if (distBg >= 0.0f) { + if (isPi) { + registry.fill(HIST("ptGenPiPerp"), mcPt, distBg, jet.detPtCorr, centrality); + } else if (isKa) { + registry.fill(HIST("ptGenKaPerp"), mcPt, distBg, jet.detPtCorr, centrality); + } else { + registry.fill(HIST("ptGenPrPerp"), mcPt, distBg, jet.detPtCorr, centrality); + } } + } - // Perp Cone - float dPhiBg1 = std::abs(mcParticle.phi() - jet.partPhiBg1); - if (dPhiBg1 > o2::constants::math::PI) - dPhiBg1 = o2::constants::math::TwoPI - dPhiBg1; - float dPhiBg2 = std::abs(mcParticle.phi() - jet.partPhiBg2); - if (dPhiBg2 > o2::constants::math::PI) - dPhiBg2 = o2::constants::math::TwoPI - dPhiBg2; - - float distBg1Sq = dEta * dEta + dPhiBg1 * dPhiBg1; - float distBg2Sq = dEta * dEta + dPhiBg2 * dPhiBg2; - - float distBgGen = -1.0f; - if (distBg1Sq < maxR2) - distBgGen = std::sqrt(distBg1Sq); - else if (distBg2Sq < maxR2) - distBgGen = std::sqrt(distBg2Sq); - - if (distBgGen >= 0.0f) { - if (isPi) - registry.fill(HIST("ptGenPiPerp"), mcPt, distBgGen, jet.detPtCorr, centrality); - else if (isKa) - registry.fill(HIST("ptGenKaPerp"), mcPt, distBgGen, jet.detPtCorr, centrality); - else if (isPr) - registry.fill(HIST("ptGenPrPerp"), mcPt, distBgGen, jet.detPtCorr, centrality); + for (const auto& cone : randomCones) { + const float dEta = mcParticle.eta() - cone.eta; + if (std::abs(dEta) >= distanceMax) { + continue; + } + const float dPhi = absDeltaPhi(mcParticle.phi(), cone.phi); + const float distSq = dEta * dEta + dPhi * dPhi; + if (distSq >= maxR2) { + continue; + } + const float distance = std::sqrt(distSq); + if (isPi) { + registry.fill(HIST("ptGenPiRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); + } else if (isKa) { + registry.fill(HIST("ptGenKaRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); + } else { + registry.fill(HIST("ptGenPrRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); } } } - // Numerator + // Reconstructed tracks for (const auto& jetTrack : jetTracks) { - if (!jetderiveddatautilities::selectTrack(jetTrack, trackSelection)) + const bool passesNominalBit = jetderiveddatautilities::selectTrack(jetTrack, trackSelection); + const bool passesWideDcaBit = jetderiveddatautilities::selectTrack(jetTrack, trackSelectionForDca); + if (!passesNominalBit && !passesWideDcaBit) { continue; + } - auto track = jetTrack.track_as(); - - if (std::abs(track.eta()) > etaTrUp) - continue; - if (track.tpcNClsCrossedRows() < nclcrossTpcMin) + const auto track = jetTrack.track_as(); + if (!passCommonTrackQuality(track)) { continue; - if (std::abs(track.dcaXY()) > dcaxyCutMax) - continue; - if (track.itsChi2NCl() > chi2ItsMax) - continue; - if (track.tpcChi2NCl() > chi2TpcMax) + } + + const bool passesNominal = passesNominalBit && passNominalTrackDca(track); + const bool passesWideDca = passesWideDcaBit; + if (!passesNominal && !passesWideDca) { continue; - if (track.tpcNClsFound() < nclTpcMin) + } + + if (!jetTrack.has_mcParticle()) { continue; - if (track.itsNCls() < nclItsMin) + } + const auto mcParticle = jetTrack.mcParticle_as(); + if (mcParticle.mcCollisionId() != reducedMcCollisionId) { continue; + } - if (!track.has_mcParticle()) + const int absPdg = std::abs(mcParticle.pdgCode()); + const bool isTruePi = absPdg == PDG_t::kPiPlus; + const bool isTrueKa = absPdg == PDG_t::kKPlus; + const bool isTruePr = absPdg == PDG_t::kProton; + if (!isTruePi && !isTrueKa && !isTruePr) { continue; - auto mcParticle = track.mcParticle(); + } - if (mcParticle.mcCollisionId() != collision.globalIndex()) - continue; + constexpr int ProducedByDecay = 4; // ROOT TMCProcess::kPDecay + const bool isPrimary = mcParticle.isPhysicalPrimary(); + const bool isDecay = !isPrimary && mcParticle.getProcess() == ProducedByDecay; + + const float recoPt = track.pt(); + const float mcPt = mcParticle.pt(); + const float dcaXy = track.dcaXY(); + const bool hasTof = track.hasTOF(); + + const bool passTpcPi = track.tpcNSigmaPi() > tpcNSigmaPiMin && + track.tpcNSigmaPi() < tpcNSigmaPiMax; + const bool passTpcPr = track.tpcNSigmaPr() > tpcNSigmaPrMin && + track.tpcNSigmaPr() < tpcNSigmaPrMax; + const bool isRecoPi = passSpeciesKinematics(recoPt, track.eta(), PionMass) && + passTpcPi && hasTof && + std::abs(track.tofNSigmaPi()) < nSigmaTofCut; + const bool isRecoPr = passSpeciesKinematics(recoPt, track.eta(), ProtonMass) && + passTpcPr && hasTof && + std::abs(track.tofNSigmaPr()) < nSigmaTofCut; + + float originBin = 2.5f; + if (isPrimary) { + originBin = 0.5f; + } else if (isDecay) { + originBin = 1.5f; + } - registry.fill(HIST("ptResolution"), track.pt(), track.pt() - mcParticle.pt()); + // DCA templates and wide->global selection transfer + auto fillInclusiveDcaAndTransfer = [&]() { + if (isRecoPi && isTruePi) { + if (passesWideDca) { + if (isPrimary) { + registry.fill(HIST("dcaPrimPiInc"), recoPt, dcaXy, centrality); + } else if (isDecay) { + registry.fill(HIST("dcaDecayPiInc"), recoPt, dcaXy, centrality); + } else { + registry.fill(HIST("dcaMatPiInc"), recoPt, dcaXy, centrality); + } + registry.fill(HIST("originCountWideInc"), recoPt, centrality, 0.5f, originBin); + } + if (passesNominal) { + registry.fill(HIST("originCountGlobalInc"), recoPt, centrality, 0.5f, originBin); + } + } - if (std::abs(mcParticle.y()) >= mcRapidityMax) - continue; + if (isRecoPr && isTruePr) { + if (passesWideDca) { + if (isPrimary) { + registry.fill(HIST("dcaPrimPrInc"), recoPt, dcaXy, centrality); + } else if (isDecay) { + registry.fill(HIST("dcaDecayPrInc"), recoPt, dcaXy, centrality); + } else { + registry.fill(HIST("dcaMatPrInc"), recoPt, dcaXy, centrality); + } + registry.fill(HIST("originCountWideInc"), recoPt, centrality, 1.5f, originBin); + } + if (passesNominal) { + registry.fill(HIST("originCountGlobalInc"), recoPt, centrality, 1.5f, originBin); + } + } + }; + fillInclusiveDcaAndTransfer(); - int absPdg = std::abs(mcParticle.pdgCode()); - bool isTruePi = (absPdg == PDG_t::kPiPlus); - bool isTrueKa = (absPdg == PDG_t::kKPlus); - bool isTruePr = (absPdg == PDG_t::kProton); - if (!isTruePi && !isTrueKa && !isTruePr) - continue; + for (const auto& jet : validJets) { + const float dEtaDet = track.eta() - jet.detEta; + + const float dPhiJetDet = absDeltaPhi(track.phi(), jet.detPhi); + const float distJetDetSq = dEtaDet * dEtaDet + dPhiJetDet * dPhiJetDet; + const bool inJetDet = distJetDetSq < maxR2; + const float distJetDet = inJetDet ? std::sqrt(distJetDetSq) : -1.0f; + + const float dPhiBg1Det = absDeltaPhi(track.phi(), jet.detPhiBg1); + const float dPhiBg2Det = absDeltaPhi(track.phi(), jet.detPhiBg2); + const float distBg1DetSq = dEtaDet * dEtaDet + dPhiBg1Det * dPhiBg1Det; + const float distBg2DetSq = dEtaDet * dEtaDet + dPhiBg2Det * dPhiBg2Det; + + float distBgDet = -1.0f; + if (distBg1DetSq < maxR2) { + distBgDet = std::sqrt(distBg1DetSq); + } else if (distBg2DetSq < maxR2) { + distBgDet = std::sqrt(distBg2DetSq); + } - const int producedByDecay = 4; - bool isPrimary = mcParticle.isPhysicalPrimary(); - bool isDecay = !isPrimary && (mcParticle.getProcess() == producedByDecay); - bool isMaterial = !isPrimary && !isDecay; + if (inJetDet) { + if (isRecoPi && isTruePi) { + if (passesWideDca) { + if (isPrimary) { + registry.fill(HIST("dcaPrimPi"), recoPt, dcaXy, distJetDet, jet.detPtCorr, centrality); + } else if (isDecay) { + registry.fill(HIST("dcaDecayPi"), recoPt, dcaXy, distJetDet, jet.detPtCorr, centrality); + } else { + registry.fill(HIST("dcaMatPi"), recoPt, dcaXy, distJetDet, jet.detPtCorr, centrality); + } + registry.fill(HIST("originCountWideJet"), recoPt, distJetDet, jet.detPtCorr, centrality, 0.5f, originBin); + } + if (passesNominal) { + registry.fill(HIST("originCountGlobalJet"), recoPt, distJetDet, jet.detPtCorr, centrality, 0.5f, originBin); + } + } + + if (isRecoPr && isTruePr) { + if (passesWideDca) { + if (isPrimary) { + registry.fill(HIST("dcaPrimPr"), recoPt, dcaXy, distJetDet, jet.detPtCorr, centrality); + } else if (isDecay) { + registry.fill(HIST("dcaDecayPr"), recoPt, dcaXy, distJetDet, jet.detPtCorr, centrality); + } else { + registry.fill(HIST("dcaMatPr"), recoPt, dcaXy, distJetDet, jet.detPtCorr, centrality); + } + registry.fill(HIST("originCountWideJet"), recoPt, distJetDet, jet.detPtCorr, centrality, 1.5f, originBin); + } + if (passesNominal) { + registry.fill(HIST("originCountGlobalJet"), recoPt, distJetDet, jet.detPtCorr, centrality, 1.5f, originBin); + } + } + } - float mcPt = mcParticle.pt(); - float recoPt = track.pt(); - float dcaXy = track.dcaXY(); - bool hasTof = track.hasTOF(); + if (distBgDet >= 0.0f) { + if (isRecoPi && isTruePi) { + if (passesWideDca) { + if (isPrimary) { + registry.fill(HIST("dcaPrimPiPerp"), recoPt, dcaXy, distBgDet, jet.detPtCorr, centrality); + } else if (isDecay) { + registry.fill(HIST("dcaDecayPiPerp"), recoPt, dcaXy, distBgDet, jet.detPtCorr, centrality); + } else { + registry.fill(HIST("dcaMatPiPerp"), recoPt, dcaXy, distBgDet, jet.detPtCorr, centrality); + } + registry.fill(HIST("originCountWidePerp"), recoPt, distBgDet, jet.detPtCorr, centrality, 0.5f, originBin); + } + if (passesNominal) { + registry.fill(HIST("originCountGlobalPerp"), recoPt, distBgDet, jet.detPtCorr, centrality, 0.5f, originBin); + } + } - // PID flag - // TPC nsigma - bool passTpcPiOnly = (track.tpcNSigmaPi() > tpcNSigmaPiMin && track.tpcNSigmaPi() < tpcNSigmaPiMax); - bool passTpcPrOnly = (track.tpcNSigmaPr() > tpcNSigmaPrMin && track.tpcNSigmaPr() < tpcNSigmaPrMax); + if (isRecoPr && isTruePr) { + if (passesWideDca) { + if (isPrimary) { + registry.fill(HIST("dcaPrimPrPerp"), recoPt, dcaXy, distBgDet, jet.detPtCorr, centrality); + } else if (isDecay) { + registry.fill(HIST("dcaDecayPrPerp"), recoPt, dcaXy, distBgDet, jet.detPtCorr, centrality); + } else { + registry.fill(HIST("dcaMatPrPerp"), recoPt, dcaXy, distBgDet, jet.detPtCorr, centrality); + } + registry.fill(HIST("originCountWidePerp"), recoPt, distBgDet, jet.detPtCorr, centrality, 1.5f, originBin); + } + if (passesNominal) { + registry.fill(HIST("originCountGlobalPerp"), recoPt, distBgDet, jet.detPtCorr, centrality, 1.5f, originBin); + } + } + } + } - // TPC + TOF - bool isRecoPi = passTpcPiOnly && hasTof && (std::abs(track.tofNSigmaPi()) < nSigmaTofCut); - bool isRecoPr = passTpcPrOnly && hasTof && (std::abs(track.tofNSigmaPr()) < nSigmaTofCut); + for (const auto& cone : randomCones) { + const float dEta = track.eta() - cone.eta; + if (std::abs(dEta) >= distanceMax) { + continue; + } + const float dPhi = absDeltaPhi(track.phi(), cone.phi); + const float distSq = dEta * dEta + dPhi * dPhi; + if (distSq >= maxR2) { + continue; + } + const float distance = std::sqrt(distSq); + + if (isRecoPi && isTruePi) { + if (passesWideDca) { + if (isPrimary) { + registry.fill(HIST("dcaPrimPiRand"), recoPt, dcaXy, distance, cone.referenceJetPtCorr, centrality); + } else if (isDecay) { + registry.fill(HIST("dcaDecayPiRand"), recoPt, dcaXy, distance, cone.referenceJetPtCorr, centrality); + } else { + registry.fill(HIST("dcaMatPiRand"), recoPt, dcaXy, distance, cone.referenceJetPtCorr, centrality); + } + registry.fill(HIST("originCountWideRand"), recoPt, distance, cone.referenceJetPtCorr, centrality, 0.5f, originBin); + } + if (passesNominal) { + registry.fill(HIST("originCountGlobalRand"), recoPt, distance, cone.referenceJetPtCorr, centrality, 0.5f, originBin); + } + } - if (isPrimary) { - // Inclusive - // TPC Efficiency - if (isTruePi) - registry.fill(HIST("effNumTpcPiInc"), mcPt, centrality); - if (isTrueKa) - registry.fill(HIST("effNumTpcKaInc"), mcPt, centrality); - if (isTruePr) - registry.fill(HIST("effNumTpcPrInc"), mcPt, centrality); - - // TOF Matching Efficiency - if (isTruePi && hasTof) + if (isRecoPr && isTruePr) { + if (passesWideDca) { + if (isPrimary) { + registry.fill(HIST("dcaPrimPrRand"), recoPt, dcaXy, distance, cone.referenceJetPtCorr, centrality); + } else if (isDecay) { + registry.fill(HIST("dcaDecayPrRand"), recoPt, dcaXy, distance, cone.referenceJetPtCorr, centrality); + } else { + registry.fill(HIST("dcaMatPrRand"), recoPt, dcaXy, distance, cone.referenceJetPtCorr, centrality); + } + registry.fill(HIST("originCountWideRand"), recoPt, distance, cone.referenceJetPtCorr, centrality, 1.5f, originBin); + } + if (passesNominal) { + registry.fill(HIST("originCountGlobalRand"), recoPt, distance, cone.referenceJetPtCorr, centrality, 1.5f, originBin); + } + } + } + + // Tracking and TOF-matching efficiency numerators + if (!passesNominal || !isPrimary) { + continue; + } + if (!passTruthKinematics(mcParticle)) { + continue; + } + + registry.fill(HIST("ptResolution"), track.pt(), track.pt() - mcPt); + + if (isTruePi) { + registry.fill(HIST("effNumTpcPiInc"), mcPt, centrality); + if (hasTof) { registry.fill(HIST("effNumTofPiInc"), mcPt, centrality); - if (isTrueKa && hasTof) + } + if (isRecoPi) { + registry.fill(HIST("effNumPidPiInc"), mcPt, centrality); + } + } else if (isTrueKa) { + registry.fill(HIST("effNumTpcKaInc"), mcPt, centrality); + if (hasTof) { registry.fill(HIST("effNumTofKaInc"), mcPt, centrality); - if (isTruePr && hasTof) + } + } else { + registry.fill(HIST("effNumTpcPrInc"), mcPt, centrality); + if (hasTof) { registry.fill(HIST("effNumTofPrInc"), mcPt, centrality); - - // Purity DCA - if (isRecoPi && isTruePi) - registry.fill(HIST("dcaPrimPiInc"), recoPt, dcaXy, centrality); - if (isRecoPr && isTruePr) - registry.fill(HIST("dcaPrimPrInc"), recoPt, dcaXy, centrality); - } else if (isDecay) { - if (isRecoPi && isTruePi) - registry.fill(HIST("dcaDecayPiInc"), recoPt, dcaXy, centrality); - if (isRecoPr && isTruePr) - registry.fill(HIST("dcaDecayPrInc"), recoPt, dcaXy, centrality); - } else if (isMaterial) { - if (isRecoPi && isTruePi) - registry.fill(HIST("dcaMatPiInc"), recoPt, dcaXy, centrality); - if (isRecoPr && isTruePr) - registry.fill(HIST("dcaMatPrInc"), recoPt, dcaXy, centrality); + } + if (isRecoPr) { + registry.fill(HIST("effNumPidPrInc"), mcPt, centrality); + } } for (const auto& jet : validJets) { - float dEta = mcParticle.eta() - jet.partEta; - - // --- Jet In-cone distance check --- - float dPhiJet = std::abs(mcParticle.phi() - jet.partPhi); - if (dPhiJet > o2::constants::math::PI) - dPhiJet = o2::constants::math::TwoPI - dPhiJet; - - float distJetReco = std::sqrt(dEta * dEta + dPhiJet * dPhiJet); - bool inJet = (distJetReco < distanceMax); - - // --- Perp Cone distance check --- - float dPhiBg1 = std::abs(mcParticle.phi() - jet.partPhiBg1); - if (dPhiBg1 > o2::constants::math::PI) - dPhiBg1 = o2::constants::math::TwoPI - dPhiBg1; - float dPhiBg2 = std::abs(mcParticle.phi() - jet.partPhiBg2); - if (dPhiBg2 > o2::constants::math::PI) - dPhiBg2 = o2::constants::math::TwoPI - dPhiBg2; - - float distBg1Sq = dEta * dEta + dPhiBg1 * dPhiBg1; - float distBg2Sq = dEta * dEta + dPhiBg2 * dPhiBg2; - - float distBgReco = -1.0f; - if (distBg1Sq < maxR2) - distBgReco = std::sqrt(distBg1Sq); - else if (distBg2Sq < maxR2) - distBgReco = std::sqrt(distBg2Sq); - - bool inPerp = (distBgReco >= 0.0f); - - if (inJet) { - if (isPrimary) { - // --- Jet In-cone 分子 --- - // TPC Efficiency - if (isTruePi) - registry.fill(HIST("effNumTpcPiJet"), mcPt, distJetReco, jet.detPtCorr, centrality); - if (isTrueKa) - registry.fill(HIST("effNumTpcKaJet"), mcPt, distJetReco, jet.detPtCorr, centrality); - if (isTruePr) - registry.fill(HIST("effNumTpcPrJet"), mcPt, distJetReco, jet.detPtCorr, centrality); - - // TOF Matching Efficiency - if (isTruePi && hasTof) - registry.fill(HIST("effNumTofPiJet"), mcPt, distJetReco, jet.detPtCorr, centrality); - if (isTrueKa && hasTof) - registry.fill(HIST("effNumTofKaJet"), mcPt, distJetReco, jet.detPtCorr, centrality); - if (isTruePr && hasTof) - registry.fill(HIST("effNumTofPrJet"), mcPt, distJetReco, jet.detPtCorr, centrality); - - // Purity DCA - if (isRecoPi && isTruePi) - registry.fill(HIST("dcaPrimPi"), recoPt, dcaXy, distJetReco, jet.detPtCorr, centrality); - if (isRecoPr && isTruePr) - registry.fill(HIST("dcaPrimPr"), recoPt, dcaXy, distJetReco, jet.detPtCorr, centrality); - } else if (isDecay) { - if (isRecoPi && isTruePi) - registry.fill(HIST("dcaDecayPi"), recoPt, dcaXy, distJetReco, jet.detPtCorr, centrality); - if (isRecoPr && isTruePr) - registry.fill(HIST("dcaDecayPr"), recoPt, dcaXy, distJetReco, jet.detPtCorr, centrality); - } else if (isMaterial) { - if (isRecoPi && isTruePi) - registry.fill(HIST("dcaMatPi"), recoPt, dcaXy, distJetReco, jet.detPtCorr, centrality); - if (isRecoPr && isTruePr) - registry.fill(HIST("dcaMatPr"), recoPt, dcaXy, distJetReco, jet.detPtCorr, centrality); + const float dEtaPart = mcParticle.eta() - jet.partEta; + + const float dPhiJetPart = absDeltaPhi(mcParticle.phi(), jet.partPhi); + const float distJetPartSq = dEtaPart * dEtaPart + dPhiJetPart * dPhiJetPart; + if (distJetPartSq < maxR2) { + const float distJetPart = std::sqrt(distJetPartSq); + if (isTruePi) { + registry.fill(HIST("effNumTpcPiJet"), mcPt, distJetPart, jet.detPtCorr, centrality); + if (hasTof) { + registry.fill(HIST("effNumTofPiJet"), mcPt, distJetPart, jet.detPtCorr, centrality); + } + if (isRecoPi) { + registry.fill(HIST("effNumPidPiJet"), mcPt, distJetPart, jet.detPtCorr, centrality); + } + } else if (isTrueKa) { + registry.fill(HIST("effNumTpcKaJet"), mcPt, distJetPart, jet.detPtCorr, centrality); + if (hasTof) { + registry.fill(HIST("effNumTofKaJet"), mcPt, distJetPart, jet.detPtCorr, centrality); + } + } else { + registry.fill(HIST("effNumTpcPrJet"), mcPt, distJetPart, jet.detPtCorr, centrality); + if (hasTof) { + registry.fill(HIST("effNumTofPrJet"), mcPt, distJetPart, jet.detPtCorr, centrality); + } + if (isRecoPr) { + registry.fill(HIST("effNumPidPrJet"), mcPt, distJetPart, jet.detPtCorr, centrality); + } + } + } + + const float dPhiBg1Part = absDeltaPhi(mcParticle.phi(), jet.partPhiBg1); + const float dPhiBg2Part = absDeltaPhi(mcParticle.phi(), jet.partPhiBg2); + const float distBg1PartSq = dEtaPart * dEtaPart + dPhiBg1Part * dPhiBg1Part; + const float distBg2PartSq = dEtaPart * dEtaPart + dPhiBg2Part * dPhiBg2Part; + + float distBgPart = -1.0f; + if (distBg1PartSq < maxR2) { + distBgPart = std::sqrt(distBg1PartSq); + } else if (distBg2PartSq < maxR2) { + distBgPart = std::sqrt(distBg2PartSq); + } + + if (distBgPart >= 0.0f) { + if (isTruePi) { + registry.fill(HIST("effNumTpcPiPerp"), mcPt, distBgPart, jet.detPtCorr, centrality); + if (hasTof) { + registry.fill(HIST("effNumTofPiPerp"), mcPt, distBgPart, jet.detPtCorr, centrality); + } + if (isRecoPi) { + registry.fill(HIST("effNumPidPiPerp"), mcPt, distBgPart, jet.detPtCorr, centrality); + } + } else if (isTrueKa) { + registry.fill(HIST("effNumTpcKaPerp"), mcPt, distBgPart, jet.detPtCorr, centrality); + if (hasTof) { + registry.fill(HIST("effNumTofKaPerp"), mcPt, distBgPart, jet.detPtCorr, centrality); + } + } else { + registry.fill(HIST("effNumTpcPrPerp"), mcPt, distBgPart, jet.detPtCorr, centrality); + if (hasTof) { + registry.fill(HIST("effNumTofPrPerp"), mcPt, distBgPart, jet.detPtCorr, centrality); + } + if (isRecoPr) { + registry.fill(HIST("effNumPidPrPerp"), mcPt, distBgPart, jet.detPtCorr, centrality); + } } } + } - if (inPerp) { - if (isPrimary) { - // Perp Cone - // TPC Efficiency - if (isTruePi) - registry.fill(HIST("effNumTpcPiPerp"), mcPt, distBgReco, jet.detPtCorr, centrality); - if (isTrueKa) - registry.fill(HIST("effNumTpcKaPerp"), mcPt, distBgReco, jet.detPtCorr, centrality); - if (isTruePr) - registry.fill(HIST("effNumTpcPrPerp"), mcPt, distBgReco, jet.detPtCorr, centrality); - - // TOF Matching Efficiency - if (isTruePi && hasTof) - registry.fill(HIST("effNumTofPiPerp"), mcPt, distBgReco, jet.detPtCorr, centrality); - if (isTrueKa && hasTof) - registry.fill(HIST("effNumTofKaPerp"), mcPt, distBgReco, jet.detPtCorr, centrality); - if (isTruePr && hasTof) - registry.fill(HIST("effNumTofPrPerp"), mcPt, distBgReco, jet.detPtCorr, centrality); - - // Purity DCA - if (isRecoPi && isTruePi) - registry.fill(HIST("dcaPrimPiPerp"), recoPt, dcaXy, distBgReco, jet.detPtCorr, centrality); - if (isRecoPr && isTruePr) - registry.fill(HIST("dcaPrimPrPerp"), recoPt, dcaXy, distBgReco, jet.detPtCorr, centrality); - } else if (isDecay) { - if (isRecoPi && isTruePi) - registry.fill(HIST("dcaDecayPiPerp"), recoPt, dcaXy, distBgReco, jet.detPtCorr, centrality); - if (isRecoPr && isTruePr) - registry.fill(HIST("dcaDecayPrPerp"), recoPt, dcaXy, distBgReco, jet.detPtCorr, centrality); - } else if (isMaterial) { - if (isRecoPi && isTruePi) - registry.fill(HIST("dcaMatPiPerp"), recoPt, dcaXy, distBgReco, jet.detPtCorr, centrality); - if (isRecoPr && isTruePr) - registry.fill(HIST("dcaMatPrPerp"), recoPt, dcaXy, distBgReco, jet.detPtCorr, centrality); + for (const auto& cone : randomCones) { + const float dEta = mcParticle.eta() - cone.eta; + if (std::abs(dEta) >= distanceMax) { + continue; + } + const float dPhi = absDeltaPhi(mcParticle.phi(), cone.phi); + const float distSq = dEta * dEta + dPhi * dPhi; + if (distSq >= maxR2) { + continue; + } + const float distance = std::sqrt(distSq); + + if (isTruePi) { + registry.fill(HIST("effNumTpcPiRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); + if (hasTof) { + registry.fill(HIST("effNumTofPiRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); + } + if (isRecoPi) { + registry.fill(HIST("effNumPidPiRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); + } + } else if (isTrueKa) { + registry.fill(HIST("effNumTpcKaRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); + if (hasTof) { + registry.fill(HIST("effNumTofKaRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); + } + } else { + registry.fill(HIST("effNumTpcPrRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); + if (hasTof) { + registry.fill(HIST("effNumTofPrRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); + } + if (isRecoPr) { + registry.fill(HIST("effNumPidPrRand"), mcPt, distance, cone.referenceJetPtCorr, centrality); } } } } } - PROCESS_SWITCH(JetShapeTask, processEfficiencyAndPurity, "process MC information for Efficiency and Purity", true); + PROCESS_SWITCH(JetShapeTask, processEfficiencyAndPurity, "process MC information for Efficiency and Purity", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; }