From 0782eb359ab467966d8e5948b3f7a61b39a9dd70 Mon Sep 17 00:00:00 2001 From: Changhwan Choi Date: Sun, 15 Nov 2026 19:41:56 +0900 Subject: [PATCH] Added weighted-MC processes for jet background analysis --- PWGJE/Core/JetFindingUtilities.h | 48 +++-- PWGJE/Tasks/jetBackgroundAnalysis.cxx | 285 ++++++++++++++++++-------- 2 files changed, 234 insertions(+), 99 deletions(-) diff --git a/PWGJE/Core/JetFindingUtilities.h b/PWGJE/Core/JetFindingUtilities.h index 9f0d76c9c08..8dd2356d666 100644 --- a/PWGJE/Core/JetFindingUtilities.h +++ b/PWGJE/Core/JetFindingUtilities.h @@ -336,6 +336,36 @@ void findJets(JetFinder& jetFinder, std::vector const& input } } +/** + * Applies the particle status and charge selections shared by jet finding and background QA. + * Kinematic acceptance cuts are applied to the input particle table separately. + */ +template +bool isParticleSelected(T const& particle, const std::string& particleSelection, int jetTypeParticleLevel, U const& pdgDatabase) +{ + if (particleSelection == "PhysicalPrimary" && !particle.isPhysicalPrimary()) { + return false; + } else if (particleSelection == "HepMCStatus" && particle.getHepMCStatusCode() != 1) { + return false; + } else if (particleSelection == "GenStatus" && particle.getGenStatusCode() != 1) { + return false; + } else if (particleSelection == "PhysicalPrimaryAndHepMCStatus" && (!particle.isPhysicalPrimary() || particle.getHepMCStatusCode() != 1)) { + return false; + } + if (std::isinf(particle.eta())) { + return false; + } + auto pdgParticle = pdgDatabase->GetParticle(particle.pdgCode()); + auto pdgCharge = pdgParticle ? std::abs(pdgParticle->Charge()) : -1.0; + if (jetTypeParticleLevel == static_cast(JetType::charged) && pdgCharge < 3.0) { + return false; + } + if (jetTypeParticleLevel == static_cast(JetType::neutral) && pdgCharge != 0.0) { + return false; + } + return true; +} + /** * Adds particles to a fastjet inputParticles list * @@ -350,26 +380,10 @@ template void analyseParticles(std::vector& inputParticles, const std::string& particleSelection, int jetTypeParticleLevel, T const& particles, o2::framework::Service pdgDatabase, const U* candidate = nullptr) { for (auto& particle : particles) { - if (particleSelection == "PhysicalPrimary" && !particle.isPhysicalPrimary()) { // CHECK : Does this exclude the HF hadron? - continue; - } else if (particleSelection == "HepMCStatus" && particle.getHepMCStatusCode() != 1) { // do we need isPhysicalPrimary as well? Note: Might give unforseen results if the generator isnt PYTHIA - continue; - } else if (particleSelection == "GenStatus" && particle.getGenStatusCode() != 1) { - continue; - } else if (particleSelection == "PhysicalPrimaryAndHepMCStatus" && (!particle.isPhysicalPrimary() || particle.getHepMCStatusCode() != 1)) { - continue; - } - if (std::isinf(particle.eta())) { + if (!isParticleSelected(particle, particleSelection, jetTypeParticleLevel, pdgDatabase)) { continue; } auto pdgParticle = pdgDatabase->GetParticle(particle.pdgCode()); - auto pdgCharge = pdgParticle ? std::abs(pdgParticle->Charge()) : -1.0; - if (jetTypeParticleLevel == static_cast(JetType::charged) && pdgCharge < 3.0) { - continue; - } - if (jetTypeParticleLevel == static_cast(JetType::neutral) && pdgCharge != 0.0) { - continue; - } if constexpr (jetcandidateutilities::isMcCandidate() && !jetv0utilities::isV0McCandidate()) { if (candidate != nullptr) { if ((*candidate).mcParticleId() == particle.globalIndex()) { diff --git a/PWGJE/Tasks/jetBackgroundAnalysis.cxx b/PWGJE/Tasks/jetBackgroundAnalysis.cxx index 6a5a8aaaae3..04b10d6ca68 100644 --- a/PWGJE/Tasks/jetBackgroundAnalysis.cxx +++ b/PWGJE/Tasks/jetBackgroundAnalysis.cxx @@ -16,6 +16,7 @@ /// \brief This file contains the implementation for the Charged Jet v2 analysis in the ALICE experiment #include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/Core/JetFindingUtilities.h" #include "PWGJE/DataModel/Jet.h" #include "PWGJE/DataModel/JetReducedData.h" #include "PWGJE/DataModel/JetSubtraction.h" @@ -30,6 +31,7 @@ #include #include #include +#include #include #include @@ -39,8 +41,6 @@ #include #include -#include - using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -56,6 +56,7 @@ struct JetBackgroundAnalysisTask { Configurable trackOccupancyInTimeRangeMax{"trackOccupancyInTimeRangeMax", 999999, "maximum track occupancy of collisions in neighbouring collisions in a given time range; only applied to reconstructed collisions (data and mcd jets), not mc collisions (mcp jets)"}; Configurable trackOccupancyInTimeRangeMin{"trackOccupancyInTimeRangeMin", -999999, "minimum track occupancy of collisions in neighbouring collisions in a given time range; only applied to reconstructed collisions (data and mcd jets), not mc collisions (mcp jets)"}; Configurable skipMBGapEvents{"skipMBGapEvents", false, "flag to choose to reject min. bias gap events"}; + Configurable applyMCDSelectionToMCP{"applyMCDSelectionToMCP", false, "require the first matched reconstructed collision to pass MCD event selection in MCP processes"}; Configurable nBinsFluct{"nBinsFluct", 1000, "number of bins for flucuations axes"}; Configurable deltaPtMinMaxFluct{"deltaPtMinMaxFluct", 200, "maximum and minimum delta pt value for fluctuation axes"}; @@ -64,10 +65,7 @@ struct JetBackgroundAnalysisTask { Configurable trackPtMin{"trackPtMin", 0.15, "minimum pT acceptance for tracks"}; Configurable trackPtMax{"trackPtMax", 100.0, "maximum pT acceptance for tracks"}; Configurable trackSelections{"trackSelections", "globalTracks", "set track selections"}; - - // // If weights are implemented for MCD bkg checks, the cut on pTHatMax of jets should be introduced - // Configurable pTHatMaxMCD{"pTHatMaxMCD", 999.0, "maximum fraction of hard scattering for jet acceptance in detector MC"}; - // Configurable pTHatExponent{"pTHatExponent", 6.0, "exponent of the event weight for the calculation of pTHat"}; + Configurable particleSelections{"particleSelections", "PhysicalPrimary", "particle status selection, as in the jet finder"}; Configurable randomConeR{"randomConeR", 0.4, "size of random Cone for estimating background fluctuations"}; Configurable randomConeLeadJetDeltaR{"randomConeLeadJetDeltaR", -99.0, "min distance between leading jet axis and random cone (RC) axis; if negative, min distance is set to automatic value of R_leadJet+R_RC "}; @@ -80,6 +78,22 @@ struct JetBackgroundAnalysisTask { void init(o2::framework::InitContext&) { + // Data/MCD and MCP use separate histograms and can run together. + // Within each level, weighted and unweighted processes share histograms + // and must not run together. Weighted processes subscribe to outlier tables. + if (doprocessRho && doprocessRhoMCDWeighted) { + LOGF(fatal, "Disable processRho when enabling processRhoMCDWeighted"); + } + if (doprocessRhoMCP && doprocessRhoMCPWeighted) { + LOGF(fatal, "Enable only one of processRhoMCP and processRhoMCPWeighted"); + } + if (doprocessBkgFluctuationsMCDWeighted && (doprocessBkgFluctuationsData || doprocessBkgFluctuationsMCD)) { + LOGF(fatal, "Disable processBkgFluctuationsData and processBkgFluctuationsMCD when enabling processBkgFluctuationsMCDWeighted"); + } + if (doprocessBkgFluctuationsMCP && doprocessBkgFluctuationsMCPWeighted) { + LOGF(fatal, "Enable only one of processBkgFluctuationsMCP and processBkgFluctuationsMCPWeighted"); + } + // selection settings initialisation eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(eventSelections); trackSelection = jetderiveddatautilities::initialiseTrackSelection(trackSelections); @@ -90,40 +104,48 @@ struct JetBackgroundAnalysisTask { // histogram definitions - if (doprocessRho) { - registry.add("h2_centrality_ntracks", "; centrality; N_{tracks};", {HistType::kTH2F, {{1100, 0., 110.0}, {10000, 0.0, 10000.0}}}); - registry.add("h2_ntracks_rho", "; N_{tracks}; #it{rho} (GeV/area);", {HistType::kTH2F, {{10000, 0.0, 10000.0}, {400, 0.0, 400.0}}}); - registry.add("h2_ntracks_rhom", "; N_{tracks}; #it{rho}_{m} (GeV/area);", {HistType::kTH2F, {{10000, 0.0, 10000.0}, {100, 0.0, 100.0}}}); - registry.add("h2_centrality_rho", "; centrality; #it{rho} (GeV/area);", {HistType::kTH2F, {{1100, 0., 110.}, {400, 0., 400.0}}}); - registry.add("h2_centrality_rhom", ";centrality; #it{rho}_{m} (GeV/area)", {HistType::kTH2F, {{1100, 0., 110.}, {100, 0., 100.0}}}); + if (doprocessRho || doprocessRhoMCDWeighted) { + const bool doSumw2 = doprocessRhoMCDWeighted; + registry.add("h2_centrality_ntracks", "; centrality; N_{tracks};", {HistType::kTH2F, {{1100, 0., 110.0}, {10000, 0.0, 10000.0}}}, doSumw2); + registry.add("h2_ntracks_rho", "; N_{tracks}; #it{rho} (GeV/area);", {HistType::kTH2F, {{10000, 0.0, 10000.0}, {400, 0.0, 400.0}}}, doSumw2); + registry.add("h2_ntracks_rhom", "; N_{tracks}; #it{rho}_{m} (GeV/area);", {HistType::kTH2F, {{10000, 0.0, 10000.0}, {100, 0.0, 100.0}}}, doSumw2); + registry.add("h2_centrality_rho", "; centrality; #it{rho} (GeV/area);", {HistType::kTH2F, {{1100, 0., 110.}, {400, 0., 400.0}}}, doSumw2); + registry.add("h2_centrality_rhom", ";centrality; #it{rho}_{m} (GeV/area)", {HistType::kTH2F, {{1100, 0., 110.}, {100, 0., 100.0}}}, doSumw2); } - if (doprocessRhoMCP) { - registry.add("h2_nparticles_rho_mcp", "; N_{tracks}; #it{rho} (GeV/area);", {HistType::kTH2F, {{10000, 0.0, 10000.0}, {400, 0.0, 400.0}}}); - registry.add("h2_nparticles_rhom_mcp", "; N_{tracks}; #it{rho}_{m} (GeV/area);", {HistType::kTH2F, {{10000, 0.0, 10000.0}, {100, 0.0, 100.0}}}); + if (doprocessRhoMCP || doprocessRhoMCPWeighted) { + const bool doSumw2 = doprocessRhoMCPWeighted; + registry.add("h2_nparticles_rho_mcp", "; N_{tracks}; #it{rho} (GeV/area);", {HistType::kTH2F, {{10000, 0.0, 10000.0}, {400, 0.0, 400.0}}}, doSumw2); + registry.add("h2_nparticles_rhom_mcp", "; N_{tracks}; #it{rho}_{m} (GeV/area);", {HistType::kTH2F, {{10000, 0.0, 10000.0}, {100, 0.0, 100.0}}}, doSumw2); } - if (doprocessBkgFluctuationsData || doprocessBkgFluctuationsMCD) { - registry.add("h2_centrality_rhorandomcone", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}); - registry.add("h2_centrality_rhorandomconerandomtrackdirection", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}); - registry.add("h2_centrality_rhorandomconewithoutleadingjet", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}); - registry.add("h2_centrality_rhorandomconerandomtrackdirectionwithoutoneleadingjets", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}); - registry.add("h2_centrality_rhorandomconerandomtrackdirectionwithouttwoleadingjets", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}); + if (doprocessBkgFluctuationsData || doprocessBkgFluctuationsMCD || doprocessBkgFluctuationsMCDWeighted) { + const bool doSumw2 = doprocessBkgFluctuationsMCDWeighted; + registry.add("h2_centrality_rhorandomcone", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}, doSumw2); + registry.add("h2_centrality_rhorandomconerandomtrackdirection", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}, doSumw2); + registry.add("h2_centrality_rhorandomconewithoutleadingjet", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}, doSumw2); + registry.add("h2_centrality_rhorandomconerandomtrackdirectionwithoutoneleadingjets", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}, doSumw2); + registry.add("h2_centrality_rhorandomconerandomtrackdirectionwithouttwoleadingjets", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, bkgFluctuationsAxis}}, doSumw2); } - if (doprocessBkgFluctuationsMCP) { - registry.add("h_rhorandomcone", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}); - registry.add("h_rhorandomconerandomtrackdirection", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}); - registry.add("h_rhorandomconewithoutleadingjet", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}); - registry.add("h_rhorandomconerandomtrackdirectionwithoutoneleadingjets", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}); - registry.add("h_rhorandomconerandomtrackdirectionwithouttwoleadingjets", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}); + if (doprocessBkgFluctuationsMCP || doprocessBkgFluctuationsMCPWeighted) { + const bool doSumw2 = doprocessBkgFluctuationsMCPWeighted; + registry.add("h_rhorandomcone", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}, doSumw2); + registry.add("h_rhorandomconerandomtrackdirection", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}, doSumw2); + registry.add("h_rhorandomconewithoutleadingjet", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}, doSumw2); + registry.add("h_rhorandomconerandomtrackdirectionwithoutoneleadingjets", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}, doSumw2); + registry.add("h_rhorandomconerandomtrackdirectionwithouttwoleadingjets", "; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH1F, {bkgFluctuationsAxis}}, doSumw2); } } Filter trackCuts = (aod::jtrack::pt >= trackPtMin && aod::jtrack::pt < trackPtMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax); - Filter particleCuts = (aod::jmcparticle::eta > trackEtaMin && aod::jmcparticle::eta < trackEtaMax); + Filter particleCuts = (aod::jmcparticle::pt >= trackPtMin && aod::jmcparticle::pt < trackPtMax && aod::jmcparticle::eta > trackEtaMin && aod::jmcparticle::eta < trackEtaMax); Filter eventCuts = (nabs(aod::jcollision::posZ) < vertexZCut); Filter mcEventCuts = (nabs(aod::jmccollision::posZ) < vertexZCut); + PresliceUnsorted collisionsPerMCPCollision = aod::jmccollisionlb::mcCollisionId; + Preslice particlesPerMCPCollision = aod::jmcparticle::mcCollisionId; + Preslice jetsPerMCPCollision = aod::jet::mcCollisionId; + template bool trackIsInJet(TTracks const& track, TJets const& jet) { @@ -136,7 +158,7 @@ struct JetBackgroundAnalysisTask { } template - void bkgFluctuationsRandomCone(TCollisions const& collision, TJets const& jets, TTracks const& tracks, float centrality) + void bkgFluctuationsRandomCone(TCollisions const& collision, TJets const& jets, TTracks const& tracks, float centrality, float eventWeight = 1.0f) { float randomConeEta = randomNumber.Uniform(trackEtaMin + randomConeR, trackEtaMax - randomConeR); float randomConePhi = randomNumber.Uniform(0.0, o2::constants::math::TwoPI); @@ -150,7 +172,7 @@ struct JetBackgroundAnalysisTask { } } } - registry.fill(HIST("h2_centrality_rhorandomcone"), centrality, randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); + registry.fill(HIST("h2_centrality_rhorandomcone"), centrality, randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho(), eventWeight); // randomised eta,phi for tracks, to assess part of fluctuations coming from statistically independently emitted particles float randomConePtRandomTrackDirection = 0; @@ -163,7 +185,7 @@ struct JetBackgroundAnalysisTask { } } } - registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirection"), centrality, randomConePtRandomTrackDirection - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); + registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirection"), centrality, randomConePtRandomTrackDirection - o2::constants::math::PI * randomConeR * randomConeR * collision.rho(), eventWeight); // removing the leading jet from the random cone const bool hasLead = jets.size() >= 1; @@ -194,7 +216,7 @@ struct JetBackgroundAnalysisTask { } } } - registry.fill(HIST("h2_centrality_rhorandomconewithoutleadingjet"), centrality, randomConePtWithoutLeadingJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); + registry.fill(HIST("h2_centrality_rhorandomconewithoutleadingjet"), centrality, randomConePtWithoutLeadingJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho(), eventWeight); // randomised eta,phi for tracks, to assess part of fluctuations coming from statistically independently emitted particles, removing tracks from 2 leading jets double randomConePtWithoutOneLeadJet = randomConePtRandomTrackDirection; @@ -219,35 +241,41 @@ struct JetBackgroundAnalysisTask { } } } - registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirectionwithoutoneleadingjets"), centrality, randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); - registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirectionwithouttwoleadingjets"), centrality, randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); + registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirectionwithoutoneleadingjets"), centrality, randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho(), eventWeight); + registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirectionwithouttwoleadingjets"), centrality, randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho(), eventWeight); } - template - void bkgFluctuationsRandomConeMCP(TCollisions const& collision, TJets const& jets, TTracks const& tracks) + template + void bkgFluctuationsRandomConeMCP(TMcCollision const& mcCollision, TJets const& jets, TTracks const& tracks, float eventWeight = 1.0f) { float randomConeEta = randomNumber.Uniform(trackEtaMin + randomConeR, trackEtaMax - randomConeR); float randomConePhi = randomNumber.Uniform(0.0, o2::constants::math::TwoPI); float randomConePt = 0; for (auto const& track : tracks) { + if (!jetfindingutilities::isParticleSelected(track, particleSelections, static_cast(JetType::charged), pdgDatabase)) { + continue; + } float dPhi = RecoDecay::constrainAngle(track.phi() - randomConePhi, (-o2::constants::math::PI)); float dEta = track.eta() - randomConeEta; if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) { randomConePt += track.pt(); } } - registry.fill(HIST("h_rhorandomcone"), randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); + registry.fill(HIST("h_rhorandomcone"), randomConePt - o2::constants::math::PI * randomConeR * randomConeR * mcCollision.rho(), eventWeight); // randomised eta,phi for tracks, to assess part of fluctuations coming from statistically independently emitted particles float randomConePtRandomTrackDirection = 0; for (auto const& track : tracks) { + if (!jetfindingutilities::isParticleSelected(track, particleSelections, static_cast(JetType::charged), pdgDatabase)) { + continue; + } float dPhi = RecoDecay::constrainAngle(randomNumber.Uniform(0.0, o2::constants::math::TwoPI) - randomConePhi, (-o2::constants::math::PI)); // ignores actual phi of track float dEta = randomNumber.Uniform(trackEtaMin, trackEtaMax) - randomConeEta; // ignores actual eta of track if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) { randomConePtRandomTrackDirection += track.pt(); } } - registry.fill(HIST("h_rhorandomconerandomtrackdirection"), randomConePtRandomTrackDirection - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); + registry.fill(HIST("h_rhorandomconerandomtrackdirection"), randomConePtRandomTrackDirection - o2::constants::math::PI * randomConeR * randomConeR * mcCollision.rho(), eventWeight); // removing the leading jet from the random cone const bool hasLead = jets.size() >= 1; @@ -268,6 +296,9 @@ struct JetBackgroundAnalysisTask { if (jetWasInCone) { randomConePtWithoutLeadingJet = 0.0; for (auto const& track : tracks) { + if (!jetfindingutilities::isParticleSelected(track, particleSelections, static_cast(JetType::charged), pdgDatabase)) { + continue; + } float dPhi = RecoDecay::constrainAngle(track.phi() - randomConePhi, (-o2::constants::math::PI)); float dEta = track.eta() - randomConeEta; if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) { @@ -276,7 +307,7 @@ struct JetBackgroundAnalysisTask { } } } - registry.fill(HIST("h_rhorandomconewithoutleadingjet"), randomConePtWithoutLeadingJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); + registry.fill(HIST("h_rhorandomconewithoutleadingjet"), randomConePtWithoutLeadingJet - o2::constants::math::PI * randomConeR * randomConeR * mcCollision.rho(), eventWeight); // randomised eta,phi for tracks, to assess part of fluctuations coming from statistically independently emitted particles, removing tracks from 2 leading jets double randomConePtWithoutOneLeadJet = randomConePtRandomTrackDirection; @@ -285,6 +316,9 @@ struct JetBackgroundAnalysisTask { randomConePtWithoutOneLeadJet = 0.0; randomConePtWithoutTwoLeadJet = 0.0; for (auto const& track : tracks) { + if (!jetfindingutilities::isParticleSelected(track, particleSelections, static_cast(JetType::charged), pdgDatabase)) { + continue; + } float dPhi = RecoDecay::constrainAngle(randomNumber.Uniform(0.0, o2::constants::math::TwoPI) - randomConePhi, (-o2::constants::math::PI)); // ignores actual phi of track float dEta = randomNumber.Uniform(trackEtaMin, trackEtaMax) - randomConeEta; // ignores actual eta of track if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) { @@ -299,93 +333,180 @@ struct JetBackgroundAnalysisTask { } } } - registry.fill(HIST("h_rhorandomconerandomtrackdirectionwithoutoneleadingjets"), randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); - registry.fill(HIST("h_rhorandomconerandomtrackdirectionwithouttwoleadingjets"), randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho()); + registry.fill(HIST("h_rhorandomconerandomtrackdirectionwithoutoneleadingjets"), randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * mcCollision.rho(), eventWeight); + registry.fill(HIST("h_rhorandomconerandomtrackdirectionwithouttwoleadingjets"), randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * mcCollision.rho(), eventWeight); } - void processRho(soa::Filtered>::iterator const& collision, soa::Filtered const& tracks) + template + bool isCollisionSelected(TCollision const& collision) { + // Also check the vertex here for the unfiltered matched-collision table. + if constexpr (isMCPMatched) { + if (!(std::abs(collision.posZ()) < vertexZCut)) { + return false; + } + } if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents)) { - return; + return false; } if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) { - return; + return false; } float centrality = checkCentFT0M ? collision.centFT0M() : collision.centFT0C(); - if (centrality < centralityMin || centralityMax < centrality) { - return; + return centrality >= centralityMin && centrality <= centralityMax; + } + + template + bool isMcCollisionSelected(TMcCollision const& mcCollision, TCollisions const& collisions) + { + if (!applyMCDSelectionToMCP) { + return true; + } + + const auto matchedCollisions = collisions.sliceBy(collisionsPerMCPCollision, mcCollision.globalIndex()); + if (matchedCollisions.size() == 0) { + return false; + } + + const auto collision = matchedCollisions.begin(); + if constexpr (isWeighted) { + if (collision.isOutlier()) { + return false; + } } + return isCollisionSelected(collision); + } + template + void fillRhoHistograms(TCollision const& collision, TTracks const& tracks, float eventWeight = 1.0f) + { + if (!isCollisionSelected(collision)) { + return; + } + float centrality = checkCentFT0M ? collision.centFT0M() : collision.centFT0C(); int nTracks = 0; for (auto const& track : tracks) { if (jetderiveddatautilities::selectTrack(track, trackSelection)) { nTracks++; } } - registry.fill(HIST("h2_centrality_ntracks"), centrality, nTracks); - registry.fill(HIST("h2_ntracks_rho"), nTracks, collision.rho()); - registry.fill(HIST("h2_ntracks_rhom"), nTracks, collision.rhoM()); - registry.fill(HIST("h2_centrality_rho"), centrality, collision.rho()); - registry.fill(HIST("h2_centrality_rhom"), centrality, collision.rhoM()); + registry.fill(HIST("h2_centrality_ntracks"), centrality, nTracks, eventWeight); + registry.fill(HIST("h2_ntracks_rho"), nTracks, collision.rho(), eventWeight); + registry.fill(HIST("h2_ntracks_rhom"), nTracks, collision.rhoM(), eventWeight); + registry.fill(HIST("h2_centrality_rho"), centrality, collision.rho(), eventWeight); + registry.fill(HIST("h2_centrality_rhom"), centrality, collision.rhoM(), eventWeight); } - PROCESS_SWITCH(JetBackgroundAnalysisTask, processRho, "QA for rho-area subtracted jets", true); - void processRhoMCP(soa::Filtered>::iterator const& mcCollision, soa::Filtered const& particles) + template + void fillRhoMCPHistograms(TMcCollision const& mcCollision, TParticles const& particles, float eventWeight = 1.0f) { - // no event selection as this is used for gen-only simulations + // Optional reconstructed-collision selection is applied in the MCP processes. int nParticles = 0; for (auto const& particle : particles) { - auto pdgParticle = pdgDatabase->GetParticle(particle.pdgCode()); - auto pdgCharge = pdgParticle ? std::abs(pdgParticle->Charge()) : -1.0; - if (pdgCharge > 0) { - nParticles++; + if (!jetfindingutilities::isParticleSelected(particle, particleSelections, static_cast(JetType::charged), pdgDatabase)) { + continue; } + nParticles++; } - registry.fill(HIST("h2_nparticles_rho_mcp"), nParticles, mcCollision.rho()); - registry.fill(HIST("h2_nparticles_rhom_mcp"), nParticles, mcCollision.rhoM()); + registry.fill(HIST("h2_nparticles_rho_mcp"), nParticles, mcCollision.rho(), eventWeight); + registry.fill(HIST("h2_nparticles_rhom_mcp"), nParticles, mcCollision.rhoM(), eventWeight); } - PROCESS_SWITCH(JetBackgroundAnalysisTask, processRhoMCP, "QA for rho-area subtracted jets at the MCP level", false); - void processBkgFluctuationsData(soa::Filtered>::iterator const& collision, soa::Join const& jets, soa::Filtered const& tracks) + template + void analyseBkgFluctuations(TCollision const& collision, TJets const& jets, TTracks const& tracks, float eventWeight = 1.0f) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents)) { - return; - } - if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) { + if (!isCollisionSelected(collision)) { return; } float centrality = checkCentFT0M ? collision.centFT0M() : collision.centFT0C(); - if (centrality < centralityMin || centralityMax < centrality) { - return; - } + bkgFluctuationsRandomCone(collision, jets, tracks, centrality, eventWeight); + } - bkgFluctuationsRandomCone(collision, jets, tracks, centrality); + void processRho(soa::Filtered>::iterator const& collision, soa::Filtered const& tracks) + { + fillRhoHistograms(collision, tracks); } - PROCESS_SWITCH(JetBackgroundAnalysisTask, processBkgFluctuationsData, "QA for random cone estimation of background fluctuations in data", false); + PROCESS_SWITCH(JetBackgroundAnalysisTask, processRho, "QA for rho-area subtracted jets", true); - void processBkgFluctuationsMCD(soa::Filtered>::iterator const& collision, soa::Join const& jets, soa::Filtered const& tracks) + void processRhoMCDWeighted(soa::Filtered>::iterator const& collision, soa::Filtered const& tracks) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents)) { + if (collision.isOutlier()) { return; } - if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) { - return; + fillRhoHistograms(collision, tracks, collision.weight()); + } + PROCESS_SWITCH(JetBackgroundAnalysisTask, processRhoMCDWeighted, "Weighted rho QA at the MCD level with outlier rejection", false); + + void processRhoMCP(soa::Filtered> const& mcCollisions, soa::Filtered const& particles, aod::JetCollisionsMCD const& collisions) + { + for (auto const& mcCollision : mcCollisions) { + if (!isMcCollisionSelected(mcCollision, collisions)) { + continue; + } + const auto particlesPerCollision = particles.sliceBy(particlesPerMCPCollision, mcCollision.globalIndex()); + fillRhoMCPHistograms(mcCollision, particlesPerCollision); } - float centrality = checkCentFT0M ? collision.centFT0M() : collision.centFT0C(); - if (centrality < centralityMin || centralityMax < centrality) { - return; + } + PROCESS_SWITCH(JetBackgroundAnalysisTask, processRhoMCP, "QA for rho-area subtracted jets at the MCP level", false); + + void processRhoMCPWeighted(soa::Filtered> const& mcCollisions, soa::Filtered const& particles, soa::Join const& collisions) + { + for (auto const& mcCollision : mcCollisions) { + if (mcCollision.isOutlier() || !isMcCollisionSelected(mcCollision, collisions)) { + continue; + } + const auto particlesPerCollision = particles.sliceBy(particlesPerMCPCollision, mcCollision.globalIndex()); + fillRhoMCPHistograms(mcCollision, particlesPerCollision, mcCollision.weight()); } + } + PROCESS_SWITCH(JetBackgroundAnalysisTask, processRhoMCPWeighted, "Weighted rho QA at the MCP level with outlier rejection", false); + + void processBkgFluctuationsData(soa::Filtered>::iterator const& collision, soa::Join const& jets, soa::Filtered const& tracks) + { + analyseBkgFluctuations(collision, jets, tracks); + } + PROCESS_SWITCH(JetBackgroundAnalysisTask, processBkgFluctuationsData, "QA for random cone estimation of background fluctuations in data", false); - bkgFluctuationsRandomCone(collision, jets, tracks, centrality); + void processBkgFluctuationsMCD(soa::Filtered>::iterator const& collision, soa::Join const& jets, soa::Filtered const& tracks) + { + analyseBkgFluctuations(collision, jets, tracks); } PROCESS_SWITCH(JetBackgroundAnalysisTask, processBkgFluctuationsMCD, "QA for random cone estimation of background fluctuations in mcd", false); - void processBkgFluctuationsMCP(soa::Filtered>::iterator const& collision, soa::Join const& jets, soa::Filtered const& tracks) + void processBkgFluctuationsMCDWeighted(soa::Filtered>::iterator const& collision, soa::Join const& jets, soa::Filtered const& tracks) + { + if (collision.isOutlier()) { + return; + } + analyseBkgFluctuations(collision, jets, tracks, collision.weight()); + } + PROCESS_SWITCH(JetBackgroundAnalysisTask, processBkgFluctuationsMCDWeighted, "Weighted random cone background fluctuations at the MCD level with outlier rejection", false); + + void processBkgFluctuationsMCP(soa::Filtered> const& mcCollisions, soa::Join const& jets, soa::Filtered const& particles, aod::JetCollisionsMCD const& collisions) { - // no event selection as this is used for gen-only simulations - bkgFluctuationsRandomConeMCP(collision, jets, tracks); + for (auto const& mcCollision : mcCollisions) { + if (!isMcCollisionSelected(mcCollision, collisions)) { + continue; + } + const auto jetsPerCollision = jets.sliceBy(jetsPerMCPCollision, mcCollision.globalIndex()); + const auto particlesPerCollision = particles.sliceBy(particlesPerMCPCollision, mcCollision.globalIndex()); + bkgFluctuationsRandomConeMCP(mcCollision, jetsPerCollision, particlesPerCollision); + } } PROCESS_SWITCH(JetBackgroundAnalysisTask, processBkgFluctuationsMCP, "QA for random cone estimation of background fluctuations in mcp", false); + + void processBkgFluctuationsMCPWeighted(soa::Filtered> const& mcCollisions, soa::Join const& jets, soa::Filtered const& particles, soa::Join const& collisions) + { + for (auto const& mcCollision : mcCollisions) { + if (mcCollision.isOutlier() || !isMcCollisionSelected(mcCollision, collisions)) { + continue; + } + const auto jetsPerCollision = jets.sliceBy(jetsPerMCPCollision, mcCollision.globalIndex()); + const auto particlesPerCollision = particles.sliceBy(particlesPerMCPCollision, mcCollision.globalIndex()); + bkgFluctuationsRandomConeMCP(mcCollision, jetsPerCollision, particlesPerCollision, mcCollision.weight()); + } + } + PROCESS_SWITCH(JetBackgroundAnalysisTask, processBkgFluctuationsMCPWeighted, "Weighted random cone background fluctuations at the MCP level with outlier rejection", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)