diff --git a/PWGJE/TableProducer/CMakeLists.txt b/PWGJE/TableProducer/CMakeLists.txt index 1c24d632022..d54513e4ebd 100644 --- a/PWGJE/TableProducer/CMakeLists.txt +++ b/PWGJE/TableProducer/CMakeLists.txt @@ -88,6 +88,11 @@ o2physics_add_dpl_workflow(jet-sv-reconstruction PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(quark-gluon-jets-producer + SOURCES quarkGluonJetsProducer.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGJECore FastJet::FastJet + COMPONENT_NAME Analysis) + endif() diff --git a/PWGJE/TableProducer/quarkGluonJetsProducer b/PWGJE/TableProducer/quarkGluonJetsProducer new file mode 100644 index 00000000000..98a44a726f3 --- /dev/null +++ b/PWGJE/TableProducer/quarkGluonJetsProducer @@ -0,0 +1,696 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file quarkGluonJetsProducer.cxx +/// \brief Produce a self-contained quark/gluon jet ML skim from JE derived data. + +/// For every accepted detector-level charged jet the task writes: +/// * one jet row, +/// * one row per reconstructed jet constituent, +/// * reconstructed PID information, +/// * MC-truth constituent information, +/// * leading and subleading quark/gluon candidates used for jet labelling. + +/// \author Aleksandra Mulewicz + +#include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/DataModel/Jet.h" +#include "PWGJE/DataModel/JetReducedData.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +namespace o2::aod +{ +namespace qgmljet +{ +DECLARE_SOA_COLUMN(EventId, eventId, uint64_t); +DECLARE_SOA_COLUMN(JetId, jetId, uint64_t); +DECLARE_SOA_COLUMN(McCollisionId, mcCollisionId, int32_t); +DECLARE_SOA_COLUMN(FlavorLabel, flavorLabel, int32_t); +DECLARE_SOA_COLUMN(LeadingPartonPdg, leadingPartonPdg, int32_t); +DECLARE_SOA_COLUMN(LeadingPartonPt, leadingPartonPt, float); +DECLARE_SOA_COLUMN(LeadingPartonDeltaR, leadingPartonDeltaR, float); +DECLARE_SOA_COLUMN(SubleadingPartonPdg, subleadingPartonPdg, int32_t); +DECLARE_SOA_COLUMN(SubleadingPartonPt, subleadingPartonPt, float); +DECLARE_SOA_COLUMN(SubleadingPartonDeltaR, subleadingPartonDeltaR, float); +DECLARE_SOA_COLUMN(NPartonsInCone, nPartonsInCone, int32_t); +DECLARE_SOA_COLUMN(LabelAmbiguous, labelAmbiguous, uint8_t); +DECLARE_SOA_COLUMN(JetPt, jetPt, float); +DECLARE_SOA_COLUMN(JetEta, jetEta, float); +DECLARE_SOA_COLUMN(JetPhi, jetPhi, float); +DECLARE_SOA_COLUMN(JetArea, jetArea, float); +DECLARE_SOA_COLUMN(JetRadius, jetRadius, float); +DECLARE_SOA_COLUMN(ZVertex, zVertex, float); +DECLARE_SOA_COLUMN(NConstituents, nConstituents, int32_t); +DECLARE_SOA_COLUMN(NSelectedConstituents, nSelectedConstituents, int32_t); +DECLARE_SOA_COLUMN(NPions, nPions, int32_t); +DECLARE_SOA_COLUMN(NKaons, nKaons, int32_t); +DECLARE_SOA_COLUMN(NProtons, nProtons, int32_t); +} + +DECLARE_SOA_TABLE(QGMLJets, "AOD", "QGMLJETS", + qgmljet::EventId, + qgmljet::JetId, + qgmljet::McCollisionId, + qgmljet::FlavorLabel, + qgmljet::LeadingPartonPdg, + qgmljet::LeadingPartonPt, + qgmljet::LeadingPartonDeltaR, + qgmljet::SubleadingPartonPdg, + qgmljet::SubleadingPartonPt, + qgmljet::SubleadingPartonDeltaR, + qgmljet::NPartonsInCone, + qgmljet::LabelAmbiguous, + qgmljet::JetPt, + qgmljet::JetEta, + qgmljet::JetPhi, + qgmljet::JetArea, + qgmljet::JetRadius, + qgmljet::ZVertex, + qgmljet::NConstituents, + qgmljet::NSelectedConstituents, + qgmljet::NPions, + qgmljet::NKaons, + qgmljet::NProtons); + +namespace qgmlconst +{ +DECLARE_SOA_COLUMN(EventId, eventId, uint64_t); +DECLARE_SOA_COLUMN(JetId, jetId, uint64_t); +DECLARE_SOA_COLUMN(ConstituentId, constituentId, uint64_t); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(P, p, float); +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(Px, px, float); +DECLARE_SOA_COLUMN(Py, py, float); +DECLARE_SOA_COLUMN(Pz, pz, float); +DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float); +DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float); +DECLARE_SOA_COLUMN(DeltaR, deltaR, float); +DECLARE_SOA_COLUMN(PtFraction, ptFraction, float); +DECLARE_SOA_COLUMN(Charge, charge, int32_t); +DECLARE_SOA_COLUMN(PassesAnalysisSelection, passesAnalysisSelection, uint8_t); +DECLARE_SOA_COLUMN(HasTOF, hasTOF, uint8_t); +DECLARE_SOA_COLUMN(DcaXY, dcaXY, float); +DECLARE_SOA_COLUMN(DcaZ, dcaZ, float); +DECLARE_SOA_COLUMN(TpcNSigmaPi, tpcNSigmaPi, float); +DECLARE_SOA_COLUMN(TpcNSigmaKa, tpcNSigmaKa, float); +DECLARE_SOA_COLUMN(TpcNSigmaPr, tpcNSigmaPr, float); +DECLARE_SOA_COLUMN(TofNSigmaPi, tofNSigmaPi, float); +DECLARE_SOA_COLUMN(TofNSigmaKa, tofNSigmaKa, float); +DECLARE_SOA_COLUMN(TofNSigmaPr, tofNSigmaPr, float); +DECLARE_SOA_COLUMN(TofBeta, tofBeta, float); +DECLARE_SOA_COLUMN(IsPion, isPion, uint8_t); +DECLARE_SOA_COLUMN(IsKaon, isKaon, uint8_t); +DECLARE_SOA_COLUMN(IsProton, isProton, uint8_t); +DECLARE_SOA_COLUMN(RecoPid, recoPid, int32_t); +DECLARE_SOA_COLUMN(TruthParticleId, truthParticleId, int64_t); +DECLARE_SOA_COLUMN(TruthPdg, truthPdg, int32_t); +DECLARE_SOA_COLUMN(TruthPt, truthPt, float); +DECLARE_SOA_COLUMN(TruthEta, truthEta, float); +DECLARE_SOA_COLUMN(TruthPhi, truthPhi, float); +DECLARE_SOA_COLUMN(IsPhysicalPrimary, isPhysicalPrimary, uint8_t); +} + +DECLARE_SOA_TABLE(QGMLConstituents, "AOD", "QGMLCONSTS", + qgmlconst::EventId, + qgmlconst::JetId, + qgmlconst::ConstituentId, + qgmlconst::Pt, + qgmlconst::P, + qgmlconst::Eta, + qgmlconst::Phi, + qgmlconst::Px, + qgmlconst::Py, + qgmlconst::Pz, + qgmlconst::DeltaEta, + qgmlconst::DeltaPhi, + qgmlconst::DeltaR, + qgmlconst::PtFraction, + qgmlconst::Charge, + qgmlconst::PassesAnalysisSelection, + qgmlconst::HasTOF, + qgmlconst::DcaXY, + qgmlconst::DcaZ, + qgmlconst::TpcNSigmaPi, + qgmlconst::TpcNSigmaKa, + qgmlconst::TpcNSigmaPr, + qgmlconst::TofNSigmaPi, + qgmlconst::TofNSigmaKa, + qgmlconst::TofNSigmaPr, + qgmlconst::TofBeta, + qgmlconst::IsPion, + qgmlconst::IsKaon, + qgmlconst::IsProton, + qgmlconst::RecoPid, + qgmlconst::TruthParticleId, + qgmlconst::TruthPdg, + qgmlconst::TruthPt, + qgmlconst::TruthEta, + qgmlconst::TruthPhi, + qgmlconst::IsPhysicalPrimary); +} + +using namespace o2; +using namespace o2::constants::math; +using namespace o2::framework; + +using MCDJetEvents = soa::Join; +using ChargedMCDJets = soa::Join; +using JetTrackRefs = soa::Join; + +using HadronTracksMC = soa::Join; + +struct QuarkGluonJetsProducer { + Produces qgMLJets; + Produces qgMLConstituents; + + Configurable eventSelections{"eventSelections", "selMCFull+NoITSROFrameBorder+IsGoodZvtxFT0vsPV", "JE derived-data event selections"}; + Configurable skipMBGapEvents{"skipMBGapEvents", true, "skip minimum-bias gap subgenerator events"}; + Configurable applyRCTSelection{"applyRCTSelection", false, "apply RCT selection (normally false for MC training)"}; + Configurable rctLabel{"rctLabel", "CBT_hadronPID", "RCT label used when applyRCTSelection=true"}; + std::vector eventSelectionBits; + + Configurable isppRefAnalysis{"isppRefAnalysis", true, "pp reference jet acceptance"}; + Configurable cfgEtaJetMax{"cfgEtaJetMax", 0.5, "maximum absolute jet eta in pp"}; + Configurable minJetPt{"minJetPt", 10.0, "minimum detector-level jet pT"}; + Configurable maxJetPt{"maxJetPt", 1000.0, "maximum detector-level jet pT"}; + Configurable rJet{"rJet", 0.4, "selected jet radius"}; + Configurable zVtx{"zVtx", 10.0, "maximum absolute collision z"}; + Configurable applyAreaCut{"applyAreaCut", false, "apply normalized jet-area cut"}; + Configurable minNormalizedJetArea{"minNormalizedJetArea", 0.6, "minimum A/(pi R^2)"}; + Configurable deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from tracking edge"}; + + Configurable requirePvContributor{"requirePvContributor", false, "require PV-contributor track"}; + Configurable minItsNclusters{"minItsNclusters", 5, "minimum ITS clusters"}; + Configurable minTpcNcrossedRows{"minTpcNcrossedRows", 70, "minimum TPC crossed rows"}; + Configurable minChiSquareTpc{"minChiSquareTpc", 0.5, "minimum TPC chi2/Ncl"}; + Configurable maxChiSquareTpc{"maxChiSquareTpc", 4.0, "maximum TPC chi2/Ncl"}; + Configurable maxChiSquareIts{"maxChiSquareIts", 36.0, "maximum ITS chi2/Ncl"}; + Configurable minPt{"minPt", 0.15, "minimum selected-track pT"}; + Configurable maxPt{"maxPt", 1000.0, "maximum selected-track pT"}; + Configurable minEta{"minEta", -0.8, "minimum selected-track eta"}; + Configurable maxEta{"maxEta", 0.8, "maximum selected-track eta"}; + Configurable maxDcaxy{"maxDcaxy", 0.2, "maximum absolute DCAxy"}; + Configurable maxDcaz{"maxDcaz", 0.1, "maximum absolute DCAz"}; + + Configurable labelAmbiguityPtFraction{"labelAmbiguityPtFraction", 0.5f, "ambiguous if opposite-flavour subleading parton exceeds this leading-pT fraction"}; + Configurable requireGeneratorPartons{"requireGeneratorPartons", true, "use only MC partons produced by the event generator for Q/G labelling"}; + Configurable partonGenStatus{"partonGenStatus", 0, "optional absolute generator-status filter for Q/G partons; 0 disables the status filter"}; + Configurable dropUnknownLabels{"dropUnknownLabels", false, "drop jets without quark/gluon label"}; + Configurable dropAmbiguousLabels{"dropAmbiguousLabels", false, "drop ambiguous labels"}; + + struct : ConfigurableGroup { + Configurable pidMethod{"pidMethod", 2, "PID method: 0 closest, 1 exclusive, 2 rejection-based"}; + Configurable rejectionSigma{"rejectionSigma", 3.0f, "minimum n-sigma distance required to reject competing PID hypotheses"}; + Configurable ptThresholdPion{"ptThresholdPion", 0.75f, "pT threshold above which pion PID requires combined TPC and TOF"}; + Configurable ptThresholdKaon{"ptThresholdKaon", 0.50f, "pT threshold above which kaon PID requires combined TPC and TOF"}; + Configurable ptThresholdProton{"ptThresholdProton", 0.75f, "pT threshold above which proton PID requires combined TPC and TOF"}; + Configurable nSigmaCut{"nSigmaCut", 2.0f, "PID acceptance cut for pion, kaon and proton hypotheses"}; + Configurable minPtPion{"minPtPion", 0.3f, "minimum pion pT"}; + Configurable maxPtPion{"maxPtPion", 4.0f, "maximum pion pT"}; + Configurable minPtKaon{"minPtKaon", 0.3f, "minimum kaon pT"}; + Configurable maxPtKaon{"maxPtKaon", 4.0f, "maximum kaon pT"}; + Configurable minPtProton{"minPtProton", 0.5f, "minimum proton pT"}; + Configurable maxPtProton{"maxPtProton", 4.0f, "maximum proton pT"}; + } cfg; + + Preslice mcParticlesPerCollision = aod::mcparticle::mcCollisionId; + + enum JetFlavor { UnknownJet = 0, + QuarkJet = 1, + GluonJet = 2 }; + enum RecoPidClass { UnidentifiedPid = 0, + PionPid = 1, + KaonPid = 2, + ProtonPid = 3 }; + + struct JetFlavorInfo { + int flavor{UnknownJet}; + int leadingPartonPdg{0}; + float leadingPartonPt{-1.0f}; + float leadingPartonDeltaR{-1.0f}; + int subleadingPartonPdg{0}; + float subleadingPartonPt{-1.0f}; + float subleadingPartonDeltaR{-1.0f}; + int nPartonsInCone{0}; + bool ambiguous{false}; + }; + + struct PidResult { + bool isPion{false}; + bool isKaon{false}; + bool isProton{false}; + }; + + void init(InitContext const&) + { + eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); + } + + template + bool hasITSLayerHit(TrackT const& track, int layer) const + { + return TESTBIT(track.itsClusterMap(), layer - 1); + } + + template + bool passedTrackSelection(TrackT const& track) + { + if (requirePvContributor && !track.isPVContributor()) { + return false; + } + if (!track.hasITS() || !track.hasTPC()) { + return false; + } + if (!hasITSLayerHit(track, 1) && !hasITSLayerHit(track, 2) && !hasITSLayerHit(track, 3)) { + return false; + } + if (track.itsNCls() < minItsNclusters || track.tpcNClsCrossedRows() < minTpcNcrossedRows) { + return false; + } + if (track.tpcChi2NCl() < minChiSquareTpc || track.tpcChi2NCl() > maxChiSquareTpc) { + return false; + } + if (track.itsChi2NCl() > maxChiSquareIts) { + return false; + } + if (track.eta() < minEta || track.eta() > maxEta || track.pt() < minPt || track.pt() > maxPt) { + return false; + } + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) { + return false; + } + return true; + } + + template + PidResult getPid(TrackT const& track) + { + constexpr int ClosestMatch = 0; + constexpr int ExclusiveMatch = 1; + constexpr int RejectionBased = 2; + constexpr double buffer = 999.0; + + const double pt = track.pt(); + + double dForPionsPi = 0.0; + double dForPionsKa = 0.0; + double dForPionsPr = 0.0; + + double dForKaonsPi = 0.0; + double dForKaonsKa = 0.0; + double dForKaonsPr = 0.0; + + double dForProtonsPi = 0.0; + double dForProtonsKa = 0.0; + double dForProtonsPr = 0.0; + + if (pt < cfg.ptThresholdPion) { + dForPionsPi = std::abs(track.tpcNSigmaPi()); + dForPionsKa = std::abs(track.tpcNSigmaKa()); + dForPionsPr = std::abs(track.tpcNSigmaPr()); + } else if (track.hasTOF()) { + dForPionsPi = std::hypot(track.tofNSigmaPi(), track.tpcNSigmaPi()); + dForPionsKa = std::hypot(track.tofNSigmaKa(), track.tpcNSigmaKa()); + dForPionsPr = std::hypot(track.tofNSigmaPr(), track.tpcNSigmaPr()); + } else { + dForPionsPi = buffer; + dForPionsKa = buffer; + dForPionsPr = buffer; + } + + if (pt < cfg.ptThresholdKaon) { + dForKaonsPi = std::abs(track.tpcNSigmaPi()); + dForKaonsKa = std::abs(track.tpcNSigmaKa()); + dForKaonsPr = std::abs(track.tpcNSigmaPr()); + } else if (track.hasTOF()) { + dForKaonsPi = std::hypot(track.tofNSigmaPi(), track.tpcNSigmaPi()); + dForKaonsKa = std::hypot(track.tofNSigmaKa(), track.tpcNSigmaKa()); + dForKaonsPr = std::hypot(track.tofNSigmaPr(), track.tpcNSigmaPr()); + } else { + dForKaonsPi = buffer; + dForKaonsKa = buffer; + dForKaonsPr = buffer; + } + + if (pt < cfg.ptThresholdProton) { + dForProtonsPi = std::abs(track.tpcNSigmaPi()); + dForProtonsKa = std::abs(track.tpcNSigmaKa()); + dForProtonsPr = std::abs(track.tpcNSigmaPr()); + } else if (track.hasTOF()) { + dForProtonsPi = std::hypot(track.tofNSigmaPi(), track.tpcNSigmaPi()); + dForProtonsKa = std::hypot(track.tofNSigmaKa(), track.tpcNSigmaKa()); + dForProtonsPr = std::hypot(track.tofNSigmaPr(), track.tpcNSigmaPr()); + } else { + dForProtonsPi = buffer; + dForProtonsKa = buffer; + dForProtonsPr = buffer; + } + + const bool isPiMatch = dForPionsPi <= cfg.nSigmaCut; + const bool isKaMatch = dForKaonsKa <= cfg.nSigmaCut; + const bool isPrMatch = dForProtonsPr <= cfg.nSigmaCut; + + PidResult result; + + if (cfg.pidMethod == ClosestMatch) { + if (isPiMatch && dForPionsPi < dForPionsKa && dForPionsPi < dForPionsPr) { + result.isPion = true; + } else if (isKaMatch && dForKaonsKa < dForKaonsPi && dForKaonsKa < dForKaonsPr) { + result.isKaon = true; + } else if (isPrMatch && dForProtonsPr < dForProtonsPi && dForProtonsPr < dForProtonsKa) { + result.isProton = true; + } + } else if (cfg.pidMethod == ExclusiveMatch) { + if (isPiMatch && !isKaMatch && !isPrMatch) { + result.isPion = true; + } else if (isKaMatch && !isPiMatch && !isPrMatch) { + result.isKaon = true; + } else if (isPrMatch && !isPiMatch && !isKaMatch) { + result.isProton = true; + } + } else if (cfg.pidMethod == RejectionBased) { + if (isPiMatch && dForPionsKa > cfg.rejectionSigma && dForPionsPr > cfg.rejectionSigma) { + result.isPion = true; + } else if (isKaMatch && dForKaonsPi > cfg.rejectionSigma && dForKaonsPr > cfg.rejectionSigma) { + result.isKaon = true; + } else if (isPrMatch && dForProtonsPi > cfg.rejectionSigma && dForProtonsKa > cfg.rejectionSigma) { + result.isProton = true; + } + } + + if (result.isPion && (pt < cfg.minPtPion || pt > cfg.maxPtPion)) { + result.isPion = false; + } + if (result.isKaon && (pt < cfg.minPtKaon || pt > cfg.maxPtKaon)) { + result.isKaon = false; + } + if (result.isProton && (pt < cfg.minPtProton || pt > cfg.maxPtProton)) { + result.isProton = false; + } + + return result; + } + + int getRecoPidClass(PidResult const& pid) const + { + if (pid.isPion) { + return PionPid; + } + if (pid.isKaon) { + return KaonPid; + } + if (pid.isProton) { + return ProtonPid; + } + return UnidentifiedPid; + } + + template + JetFlavorInfo getJetFlavorTag(JetT const& jet, ParticleRangeT const& particles) + { + JetFlavorInfo result; + const double radius = static_cast(jet.r()) / 100.0; + + for (auto const& particle : particles) { + if (requireGeneratorPartons && !particle.producedByGenerator()) { + continue; + } + if (partonGenStatus > 0 && std::abs(particle.getGenStatusCode()) != partonGenStatus) { + continue; + } + + const int absPdg = std::abs(particle.pdgCode()); + const bool isQuark = absPdg >= 1 && absPdg <= 6; + const bool isGluon = absPdg == 21; + if (!isQuark && !isGluon) { + continue; + } + + const double dEta = particle.eta() - jet.eta(); + const double dPhi = RecoDecay::constrainAngle(particle.phi() - jet.phi(), -PI); + const double dR = std::hypot(dEta, dPhi); + if (dR >= radius) { + continue; + } + + ++result.nPartonsInCone; + const float pt = particle.pt(); + if (pt > result.leadingPartonPt) { + result.subleadingPartonPdg = result.leadingPartonPdg; + result.subleadingPartonPt = result.leadingPartonPt; + result.subleadingPartonDeltaR = result.leadingPartonDeltaR; + result.leadingPartonPdg = particle.pdgCode(); + result.leadingPartonPt = pt; + result.leadingPartonDeltaR = dR; + } else if (pt > result.subleadingPartonPt) { + result.subleadingPartonPdg = particle.pdgCode(); + result.subleadingPartonPt = pt; + result.subleadingPartonDeltaR = dR; + } + } + + const int leadingAbs = std::abs(result.leadingPartonPdg); + if (leadingAbs >= 1 && leadingAbs <= 6) { + result.flavor = QuarkJet; + } else if (leadingAbs == 21) { + result.flavor = GluonJet; + } + + if (result.leadingPartonPt > 0.f && result.subleadingPartonPt > 0.f) { + const int subAbs = std::abs(result.subleadingPartonPdg); + const bool leadingQuark = leadingAbs >= 1 && leadingAbs <= 6; + const bool subQuark = subAbs >= 1 && subAbs <= 6; + result.ambiguous = (leadingQuark != subQuark) && + result.subleadingPartonPt >= labelAmbiguityPtFraction * result.leadingPartonPt; + } + return result; + } + + int getJetMcCollisionId(ChargedMCDJets::iterator const& jet) + { + for (auto const& jtrack : jet.tracks_as()) { + auto track = jtrack.track_as(); + if (!track.has_mcParticle()) { + continue; + } + auto mcParticle = track.mcParticle_as(); + return mcParticle.mcCollisionId(); + } + return -1; + } + + void process(MCDJetEvents::iterator const& collision, + aod::JetMcCollisions const&, + ChargedMCDJets const& detectorLevelJets, + JetTrackRefs const&, + HadronTracksMC const&, + aod::McParticles const& mcParticles) + { + if (!jetderiveddatautilities::selectCollision( + collision, + eventSelectionBits, + static_cast(skipMBGapEvents), + static_cast(applyRCTSelection), + static_cast(rctLabel))) { + return; + } + if (std::abs(collision.posZ()) > zVtx) { + return; + } + + for (auto const& jet : detectorLevelJets) { + const double jetRadius = static_cast(jet.r()) / 100.0; + if (std::abs(jetRadius - static_cast(rJet)) > 1.e-3) { + continue; + } + + if (isppRefAnalysis && std::abs(jet.eta()) > cfgEtaJetMax) { + continue; + } + if (!isppRefAnalysis && (std::abs(jet.eta()) + jetRadius > maxEta - deltaEtaEdge)) { + continue; + } + if (jet.pt() < minJetPt || jet.pt() > maxJetPt) { + continue; + } + const double normalizedArea = jet.area() / (PI * jetRadius * jetRadius); + if (applyAreaCut && normalizedArea < minNormalizedJetArea) { + continue; + } + + const int jetMcCollisionId = getJetMcCollisionId(jet); + JetFlavorInfo flavor; + if (jetMcCollisionId >= 0) { + auto collisionParticles = mcParticles.sliceBy(mcParticlesPerCollision, jetMcCollisionId); + flavor = getJetFlavorTag(jet, collisionParticles); + } + if (dropUnknownLabels && flavor.flavor == UnknownJet) { + continue; + } + if (dropAmbiguousLabels && flavor.ambiguous) { + continue; + } + + const uint64_t eventId = static_cast(collision.globalIndex()); + const uint64_t jetId = static_cast(jet.globalIndex()); + + int nSelected = 0; + int nPions = 0; + int nKaons = 0; + int nProtons = 0; + + for (auto const& jtrack : jet.tracks_as()) { + auto track = jtrack.track_as(); + + const bool passes = passedTrackSelection(track); + const double dEta = track.eta() - jet.eta(); + const double dPhi = RecoDecay::constrainAngle(track.phi() - jet.phi(), -PI); + const double dR = std::hypot(dEta, dPhi); + const PidResult pid = getPid(track); + const int recoPid = getRecoPidClass(pid); + + int64_t truthParticleId = -1; + int32_t truthPdg = 0; + float truthPt = std::numeric_limits::quiet_NaN(); + float truthEta = std::numeric_limits::quiet_NaN(); + float truthPhi = std::numeric_limits::quiet_NaN(); + uint8_t isPhysicalPrimary = 0; + + if (track.has_mcParticle()) { + auto truth = track.mcParticle_as(); + truthParticleId = static_cast(truth.globalIndex()); + truthPdg = truth.pdgCode(); + truthPt = truth.pt(); + truthEta = truth.eta(); + truthPhi = truth.phi(); + isPhysicalPrimary = static_cast(truth.isPhysicalPrimary()); + } + + const float missingTOF = std::numeric_limits::quiet_NaN(); + qgMLConstituents( + eventId, + jetId, + static_cast(track.globalIndex()), + static_cast(track.pt()), + static_cast(track.p()), + static_cast(track.eta()), + static_cast(track.phi()), + static_cast(track.px()), + static_cast(track.py()), + static_cast(track.pz()), + static_cast(dEta), + static_cast(dPhi), + static_cast(dR), + jet.pt() > 0.f ? static_cast(track.pt() / jet.pt()) : 0.f, + track.sign(), + static_cast(passes), + static_cast(track.hasTOF()), + static_cast(track.dcaXY()), + static_cast(track.dcaZ()), + static_cast(track.tpcNSigmaPi()), + static_cast(track.tpcNSigmaKa()), + static_cast(track.tpcNSigmaPr()), + track.hasTOF() ? static_cast(track.tofNSigmaPi()) : missingTOF, + track.hasTOF() ? static_cast(track.tofNSigmaKa()) : missingTOF, + track.hasTOF() ? static_cast(track.tofNSigmaPr()) : missingTOF, + track.hasTOF() ? static_cast(track.beta()) : missingTOF, + static_cast(pid.isPion), + static_cast(pid.isKaon), + static_cast(pid.isProton), + recoPid, + truthParticleId, + truthPdg, + truthPt, + truthEta, + truthPhi, + isPhysicalPrimary); + + if (!passes) { + continue; + } + ++nSelected; + if (pid.isPion) { + ++nPions; + } else if (pid.isKaon) { + ++nKaons; + } else if (pid.isProton) { + ++nProtons; + } + } + + qgMLJets( + eventId, + jetId, + jetMcCollisionId, + flavor.flavor, + flavor.leadingPartonPdg, + flavor.leadingPartonPt, + flavor.leadingPartonDeltaR, + flavor.subleadingPartonPdg, + flavor.subleadingPartonPt, + flavor.subleadingPartonDeltaR, + flavor.nPartonsInCone, + static_cast(flavor.ambiguous), + static_cast(jet.pt()), + static_cast(jet.eta()), + static_cast(jet.phi()), + static_cast(jet.area()), + static_cast(jetRadius), + static_cast(collision.posZ()), + static_cast(jet.tracksIds().size()), + nSelected, + nPions, + nKaons, + nProtons); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"quark-gluon-jets-producer"})}; +}