diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index 55ec23b2c44..f20bd617e33 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -403,4 +403,8 @@ if(FastJet_FOUND) SOURCES bjetCentMult.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(jet-hadrons-pid + SOURCES jetHadronsPid.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGJECore FastJet::FastJet FastJet::Contrib O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) endif() diff --git a/PWGJE/Tasks/jetHadronsPid.cxx b/PWGJE/Tasks/jetHadronsPid.cxx new file mode 100644 index 00000000000..3f3917b6575 --- /dev/null +++ b/PWGJE/Tasks/jetHadronsPid.cxx @@ -0,0 +1,660 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file jetHadronsPID.cxx +/// \brief Analysis of hadrons in jets +/// \author Leonard Lorenc, WUT Warsaw, leonard.lorenc@cern.ch +/// \author Aleksandra Mulewicz, WUT Warsaw, aleksandra.mulewicz@cern.ch +/// \author Hubert Zalewski, WUT Warsaw, hubert.zalewski@cern.ch +/// \author Janik Małgorzata malgorzata.janik@cern.ch +/// \author Daniela Ruggiano daniela.ruggiano@cern.ch + +#include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/Core/JetUtilities.h" +#include "PWGJE/DataModel/Jet.h" +#include "PWGJE/DataModel/JetReducedData.h" +#include "PWGJE/DataModel/JetSubtraction.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/DCA.h" +#include "ReconstructionDataFormats/PID.h" +#include "ReconstructionDataFormats/Track.h" + +#include + +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod; +using namespace o2::constants::physics; +using namespace o2::constants::math; +using namespace o2::framework::expressions; + +using StandardEvents = soa::Join; + +using JetEvents = soa::Join; + +using HadronTracks = soa::Join; + +using HadronTracksMC = soa::Join; + +struct jetHadronsPid { + + HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + Configurable eventSelections{"eventSelections", "sel8,vtxITSTPC,noITSROF,noTF,noPileup,goodZvtx", "choose event selection"}; + std::vector eventSelectionBits; + + Configurable isppRefAnalysis{"isppRefAnalysis", true, "Is ppRef analysis"}; + Configurable cfgEtaJetMax{"cfgEtaJetMax", 0.5, "max jet eta"}; + + Configurable rejectITSROFBorder{"rejectITSROFBorder", true, "Reject events near the ITS ROF border"}; + Configurable rejectTFBorder{"rejectTFBorder", true, "Reject events near the TF border"}; + Configurable requireVtxITSTPC{"requireVtxITSTPC", true, "Require at least one ITS-TPC matched track"}; + Configurable rejectSameBunchPileup{"rejectSameBunchPileup", true, "Reject events with same-bunch pileup collisions"}; + Configurable requireIsGoodZvtxFT0VsPV{"requireIsGoodZvtxFT0VsPV", true, "Require consistent FT0 vs PV z-vertex"}; + Configurable requireIsVertexTOFmatched{"requireIsVertexTOFmatched", false, "Require at least one vertex track matched to TOF"}; + + Configurable minJetPt{"minJetPt", 4.0, "Minimum pt of the jet after bkg subtraction"}; + Configurable maxJetPt{"maxJetPt", 1e+06, "Maximum pt of the jet after bkg subtraction"}; + Configurable rJet{"rJet", 0.4, "Jet resolution parameter R"}; + Configurable zVtx{"zVtx", 10.0, "Maximum zVertex"}; + Configurable applyAreaCut{"applyAreaCut", false, "apply area cut"}; + Configurable minNormalizedJetArea{"minNormalizedJetArea", 0.6, "Minimum normalized area cut to reject fake jets"}; + Configurable deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"}; + + Configurable requirePvContributor{"requirePvContributor", false, "require that the track is a PV contributor"}; + Configurable minItsNclusters{"minItsNclusters", 5, "minimum number of ITS clusters"}; + Configurable minTpcNcrossedRows{"minTpcNcrossedRows", 70, "minimum number of TPC crossed pad rows"}; + Configurable minChiSquareTpc{"minChiSquareTpc", 0.0, "minimum TPC chi^2/Ncls"}; + Configurable maxChiSquareTpc{"maxChiSquareTpc", 4.0, "maximum TPC chi^2/Ncls"}; + Configurable maxChiSquareIts{"maxChiSquareIts", 36.0, "maximum ITS chi^2/Ncls"}; + Configurable minPt{"minPt", 0.3, "minimum pt of the tracks"}; + Configurable maxPt{"maxPt", 4.0, "maximum pt of the tracks for PID analysis"}; + Configurable minEta{"minEta", -0.8, "minimum eta"}; + Configurable maxEta{"maxEta", +0.8, "maximum eta"}; + Configurable maxDcaxy{"maxDcaxy", 0.2, "Maximum DCAxy"}; + Configurable maxDcaz{"maxDcaz", 0.1, "Maximum DCAz"}; + Configurable maxNSigmaPid{"maxNSigmaPid", 3.0, "Maximum nSigma for TPC and TOF PID"}; + + Configurable setMCDefaultItsParams{"setMCDefaultItsParams", true, "set MC default parameters"}; + + o2::aod::ITSResponse itsResponse; + + Preslice tracksPerCollision = aod::track::collisionId; + + void init(InitContext const&) + { + if (setMCDefaultItsParams) { + itsResponse.setMCDefaultParameters(); + } + + eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); + + registryData.add("n_events", "Event counter", HistType::kTH1F, {{1, 0.5, 1.5, "N_{events}"}}); + registryData.add("n_events_raw", "All events", HistType::kTH1F, {{1, 0.5, 1.5, ""}}); + + registryData.add("jet_pt", "Jet pT ", HistType::kTH1F, {{100, 0.0, 20.0, "#it{p}_{T}^{raw} (GeV/#it{c})"}}); + registryData.add("jet_pt_subtracted", "Jet pT subtracted", HistType::kTH1F, {{200, 0.0, 200.0, "#it{p}_{T}^{sub} (GeV/#it{c})"}}); + registryData.add("jet_pt_raw_vs_sub", "Raw vs sub jet pT", HistType::kTH2F, {{200, 0, 200}, {200, 0, 200}}); + registryData.add("jet_eta", "Jet eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta_{jet}"}}); + registryData.add("jet_phi", "Jet phi", HistType::kTH1F, {{100, 0.0, TwoPI, "#phi_{jet}"}}); + registryData.add("jet_area", "Jet area", HistType::kTH1F, {{100, 0.0, 1.5, "Area"}}); + registryData.add("jet_n_constituents", "Jet multiplicity", HistType::kTH1I, {{100, 0, 100, "N_{constituents}"}}); + + registryData.add("pion_pure_tpc", "TPC Pion PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("pion_pure_tof", "TOF Pion PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("pion_pure_pt", "Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("pion_pure_eta", "Pion Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("pion_pure_dcaz", "Pion DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + registryData.add("pion_pure_dcaxy", "Pion DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + + registryData.add("kaon_pure_tpc", "TPC Kaon PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("kaon_pure_tof", "TOF Kaon PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("kaon_pure_pt", "Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("kaon_pure_eta", "Kaon Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("kaon_pure_dcaz", "Kaon DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + registryData.add("kaon_pure_dcaxy", "Kaon DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + + registryData.add("proton_pure_tpc", "TPC Proton PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("proton_pure_tof", "TOF Proton PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("proton_pure_pt", "Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("proton_pure_eta", "Proton Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("proton_pure_dcaz", "Proton DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + registryData.add("proton_pure_dcaxy", "Proton DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + + registryData.add("z_vtx", "Z-Vertex Distribution", HistType::kTH1F, {{200, -20.0, 20.0, "Z-Vertex (cm)"}}); + + registryData.add("pion_jet_tpc", "TPC Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("pion_jet_tof", "TOF Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("kaon_jet_tpc", "TPC Kaon PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("kaon_jet_tof", "TOF Kaon PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("proton_jet_tpc", "TPC Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("proton_jet_tof", "TOF Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + + registryData.add("pion_jet_pt", "Pion pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("pion_jet_eta", "Pion Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("pion_jet_dcaxy", "Pion DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("pion_jet_dcaz", "Pion DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + + registryData.add("kaon_jet_pt", "Kaon pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("kaon_jet_eta", "Kaon Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("kaon_jet_dcaxy", "Kaon DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("kaon_jet_dcaz", "Kaon DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + + registryData.add("proton_jet_pt", "Proton pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("proton_jet_eta", "Proton Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("proton_jet_dcaxy", "Proton DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("proton_jet_dcaz", "Proton DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + + registryData.add("pion_ue_tpc", "TPC Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("pion_ue_tof", "TOF Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("kaon_ue_tpc", "TPC Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("kaon_ue_tof", "TOF Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("proton_ue_tpc", "TPC Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("proton_ue_tof", "TOF Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + + registryData.add("pion_ue_pt", "Pion pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("pion_ue_eta", "Pion Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("pion_ue_dcaxy", "Pion DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("pion_ue_dcaz", "Pion DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + + registryData.add("kaon_ue_pt", "Kaon pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("kaon_ue_eta", "Kaon Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("kaon_ue_dcaxy", "Kaon DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("kaon_ue_dcaz", "Kaon DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + + registryData.add("proton_ue_pt", "Proton pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("proton_ue_eta", "Proton Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("proton_ue_dcaxy", "Proton DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("proton_ue_dcaz", "Proton DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + + registryData.add("tracks_n_in_ue", "Number of tracks in UE", HistType::kTH1I, {{100, 0, 100, "N_{tracks}"}}); + + registryData.add("rec_pion_all", "All Tracks PID'd as Pions", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_rec_pion_pt", "True Primary Pions (Reconstructed)", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_sec_pion_pt", "Secondary Pions (Reconstructed)", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("contamination_matrix_pion", "Pion PID Contamination", HistType::kTH2F, {{4000, -0.5, 3999.5, "Absolute PDG Code"}, {120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + + registryData.add("rec_kaon_all", "All Tracks PID'd as Kaons", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_rec_kaon_pt", "True Primary Kaons (Reconstructed)", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_sec_kaon_pt", "Secondary Kaons (Reconstructed)", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("contamination_matrix_kaon", "Kaon PID Contamination", HistType::kTH2F, {{4000, -0.5, 3999.5, "Absolute PDG Code"}, {120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + + registryData.add("rec_proton_all", "All Tracks PID'd as Protons", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_rec_proton_pt", "True Primary Protons (Reconstructed)", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_sec_proton_pt", "Secondary Protons (Reconstructed)", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("contamination_matrix_proton", "Proton PID Contamination", HistType::kTH2F, {{4000, -0.5, 3999.5, "Absolute PDG Code"}, {120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + + registryData.add("mc_gen_pion_pt", "Generated Primary Pions (Truth)", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_gen_kaon_pt", "Generated Primary Kaons (Truth)", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_gen_proton_pt", "Generated Primary Protons (Truth)", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + } + + void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2) + { + if (p.Mag2() < 1e-9) { + u1.SetXYZ(0, 0, 0); + u2.SetXYZ(0, 0, 0); + return; + } + u1 = p.Orthogonal(); + u2 = p.Cross(u1); + } + + template + bool hasITSLayerHit(const TrackIts& track, int layer) + { + int ibit = layer - 1; + return (track.itsClusterMap() & (1 << ibit)) != 0; + } + + template + bool passedTrackSelection(const PionTrack& 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) + return false; + if (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) + return false; + if (track.pt() < minPt || track.pt() > maxPt) + return false; + return true; + } + + void processPureTracks(StandardEvents::iterator const& collision, HadronTracks const& globalTracks) + { + if (!collision.sel8() || std::abs(collision.posZ()) > zVtx) + return; + if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) + return; + if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) + return; + if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) + return; + if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) + return; + if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) + return; + if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) + return; + + for (auto const& track : globalTracks) { + if (!passedTrackSelection(track)) + continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) + continue; + + double pt = track.pt(); + double eta = track.eta(); + double dcaxy = track.dcaXY(); + double dcaz = track.dcaZ(); + + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid); + + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid); + + const double ptThreshold = 0.8; + + if (passTpcPi) { + registryData.fill(HIST("pion_pure_tpc"), pt, track.tpcNSigmaPi()); + if (passTofPi) + registryData.fill(HIST("pion_pure_tof"), pt, track.tofNSigmaPi()); + if (pt < ptThreshold || passTofPi) { + registryData.fill(HIST("pion_pure_pt"), pt); + registryData.fill(HIST("pion_pure_eta"), eta); + registryData.fill(HIST("pion_pure_dcaxy"), dcaxy); + registryData.fill(HIST("pion_pure_dcaz"), dcaz); + } + } + + if (passTpcKa) { + registryData.fill(HIST("kaon_pure_tpc"), pt, track.tpcNSigmaKa()); + if (passTofKa) + registryData.fill(HIST("kaon_pure_tof"), pt, track.tofNSigmaKa()); + if (pt < ptThreshold || passTofKa) { + registryData.fill(HIST("kaon_pure_pt"), pt); + registryData.fill(HIST("kaon_pure_eta"), eta); + registryData.fill(HIST("kaon_pure_dcaxy"), dcaxy); + registryData.fill(HIST("kaon_pure_dcaz"), dcaz); + } + } + + if (passTpcPr) { + registryData.fill(HIST("proton_pure_tpc"), pt, track.tpcNSigmaPr()); + if (passTofPr) + registryData.fill(HIST("proton_pure_tof"), pt, track.tofNSigmaPr()); + if (pt < ptThreshold || passTofPr) { + registryData.fill(HIST("proton_pure_pt"), pt); + registryData.fill(HIST("proton_pure_eta"), eta); + registryData.fill(HIST("proton_pure_dcaxy"), dcaxy); + registryData.fill(HIST("proton_pure_dcaz"), dcaz); + } + } + } + } + PROCESS_SWITCH(jetHadronsPid, processPureTracks, "Pure Tracks Analysis", true); + + void processJets(JetEvents::iterator const& collision, + soa::Join const& jets, + soa::Join const&, + HadronTracks const& globalTracks) + { + registryData.fill(HIST("n_events_raw"), 1); + + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + return; + } + + float zVertex = collision.posZ(); + if (std::abs(zVertex) > zVtx) + return; + + registryData.fill(HIST("n_events"), 1); + registryData.fill(HIST("z_vtx"), zVertex); + + double centralRho = collision.rho(); + + int baseCollId = collision.collisionId(); + + auto collTracks = globalTracks.sliceBy(tracksPerCollision, baseCollId); + + for (auto& jet : jets) { + + if (!isppRefAnalysis && ((std::abs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) + continue; + if (isppRefAnalysis && std::abs(jet.eta()) > cfgEtaJetMax) + continue; + + double ptSub = jet.pt() - (centralRho * jet.area()); + if (ptSub < 0) + ptSub = 0.0; + + registryData.fill(HIST("jet_pt_subtracted"), ptSub); + registryData.fill(HIST("jet_pt_raw_vs_sub"), jet.pt(), ptSub); + + if (isppRefAnalysis && (jet.pt() < minJetPt || jet.pt() > maxJetPt)) + continue; + if (!isppRefAnalysis && (ptSub < minJetPt || ptSub > maxJetPt)) + continue; + + double normalizedJetArea = jet.area() / (PI * rJet * rJet); + if (applyAreaCut && normalizedJetArea < minNormalizedJetArea) + continue; + + registryData.fill(HIST("jet_eta"), jet.eta()); + registryData.fill(HIST("jet_phi"), jet.phi()); + registryData.fill(HIST("jet_area"), jet.area()); + registryData.fill(HIST("jet_pt"), jet.pt()); + + TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); + TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0); + getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2); + + int constituentCount = 0; + std::set tracksInJetsSet; + + for (auto& jtrack : jet.tracks_as>()) { + constituentCount++; + + auto track = jtrack.track_as(); + tracksInJetsSet.insert(track.index()); + + if (!passedTrackSelection(track)) + continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) + continue; + + double pt = track.pt(); + double eta = track.eta(); + double dcaxy = track.dcaXY(); + double dcaz = track.dcaZ(); + + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid); + + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid); + + const double ptThreshold = 0.8; + + if (passTpcPi) { + registryData.fill(HIST("pion_jet_tpc"), pt, track.tpcNSigmaPi()); + if (passTofPi) + registryData.fill(HIST("pion_jet_tof"), pt, track.tofNSigmaPi()); + if (pt < ptThreshold || passTofPi) { + registryData.fill(HIST("pion_jet_pt"), pt); + registryData.fill(HIST("pion_jet_eta"), eta); + registryData.fill(HIST("pion_jet_dcaxy"), dcaxy); + registryData.fill(HIST("pion_jet_dcaz"), dcaz); + } + } + + if (passTpcKa) { + registryData.fill(HIST("kaon_jet_tpc"), pt, track.tpcNSigmaKa()); + if (passTofKa) + registryData.fill(HIST("kaon_jet_tof"), pt, track.tofNSigmaKa()); + if (pt < ptThreshold || passTofKa) { + registryData.fill(HIST("kaon_jet_pt"), pt); + registryData.fill(HIST("kaon_jet_eta"), eta); + registryData.fill(HIST("kaon_jet_dcaxy"), dcaxy); + registryData.fill(HIST("kaon_jet_dcaz"), dcaz); + } + } + + if (passTpcPr) { + registryData.fill(HIST("proton_jet_tpc"), pt, track.tpcNSigmaPr()); + if (passTofPr) + registryData.fill(HIST("proton_jet_tof"), pt, track.tofNSigmaPr()); + if (pt < ptThreshold || passTofPr) { + registryData.fill(HIST("proton_jet_pt"), pt); + registryData.fill(HIST("proton_jet_eta"), eta); + registryData.fill(HIST("proton_jet_dcaxy"), dcaxy); + registryData.fill(HIST("proton_jet_dcaz"), dcaz); + } + } + } + registryData.fill(HIST("jet_n_constituents"), constituentCount); + + if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) + continue; + + int nTracksOut = 0; + + for (auto const& track : collTracks) { + + if (tracksInJetsSet.find(track.index()) == tracksInJetsSet.end()) { + if (passedTrackSelection(track)) + nTracksOut++; + } + + if (!passedTrackSelection(track)) + continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) + continue; + + double deltaEtaUe1 = track.eta() - ueAxis1.Eta(); + double deltaPhiUe1 = std::abs(RecoDecay::constrainAngle(track.phi() - ueAxis1.Phi())); + double deltaRUe1 = std::sqrt(deltaEtaUe1 * deltaEtaUe1 + deltaPhiUe1 * deltaPhiUe1); + + double deltaEtaUe2 = track.eta() - ueAxis2.Eta(); + double deltaPhiUe2 = std::abs(RecoDecay::constrainAngle(track.phi() - ueAxis2.Phi())); + double deltaRUe2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2); + + if (deltaRUe1 > rJet && deltaRUe2 > rJet) + continue; + + double pt = track.pt(); + double eta = track.eta(); + double dcaxy = track.dcaXY(); + double dcaz = track.dcaZ(); + + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid); + + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid); + + const double ptThreshold = 0.8; + + if (passTpcPi) { + registryData.fill(HIST("pion_ue_tpc"), pt, track.tpcNSigmaPi()); + if (passTofPi) + registryData.fill(HIST("pion_ue_tof"), pt, track.tofNSigmaPi()); + if (pt < ptThreshold || passTofPi) { + registryData.fill(HIST("pion_ue_pt"), pt); + registryData.fill(HIST("pion_ue_eta"), eta); + registryData.fill(HIST("pion_ue_dcaxy"), dcaxy); + registryData.fill(HIST("pion_ue_dcaz"), dcaz); + } + } + + if (passTpcKa) { + registryData.fill(HIST("kaon_ue_tpc"), pt, track.tpcNSigmaKa()); + if (passTofKa) + registryData.fill(HIST("kaon_ue_tof"), pt, track.tofNSigmaKa()); + if (pt < ptThreshold || passTofKa) { + registryData.fill(HIST("kaon_ue_pt"), pt); + registryData.fill(HIST("kaon_ue_eta"), eta); + registryData.fill(HIST("kaon_ue_dcaxy"), dcaxy); + registryData.fill(HIST("kaon_ue_dcaz"), dcaz); + } + } + + if (passTpcPr) { + registryData.fill(HIST("proton_ue_tpc"), pt, track.tpcNSigmaPr()); + if (passTofPr) + registryData.fill(HIST("proton_ue_tof"), pt, track.tofNSigmaPr()); + if (pt < ptThreshold || passTofPr) { + registryData.fill(HIST("proton_ue_pt"), pt); + registryData.fill(HIST("proton_ue_eta"), eta); + registryData.fill(HIST("proton_ue_dcaxy"), dcaxy); + registryData.fill(HIST("proton_ue_dcaz"), dcaz); + } + } + } + registryData.fill(HIST("tracks_n_in_ue"), nTracksOut); + } + } + PROCESS_SWITCH(jetHadronsPid, processJets, "Jets Analysis", true); + + void processMC(StandardEvents::iterator const& collision, HadronTracksMC const& tracks, aod::McParticles const&) + { + if (!collision.sel8() || std::abs(collision.posZ()) > zVtx) + return; + if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) + return; + if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) + return; + if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) + return; + if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) + return; + if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) + return; + if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) + return; + + const double ptThreshold = 0.8; + + for (auto const& track : tracks) { + + if (!passedTrackSelection(track)) + continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) + continue; + + double pt = track.pt(); + + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid); + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid); + + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid); + + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid); + + if (!track.has_mcParticle()) + continue; + auto const& trueParticle = track.mcParticle(); + + int pdg = std::abs(trueParticle.pdgCode()); + bool isPrimary = trueParticle.isPhysicalPrimary(); + + if (passTpcPi && (pt < ptThreshold || passTofPi)) { + registryData.fill(HIST("rec_pion_all"), pt); + registryData.fill(HIST("contamination_matrix_pion"), pdg, pt); + if (isPrimary) { + if (pdg == 211) + registryData.fill(HIST("mc_rec_pion_pt"), pt); + } else { + registryData.fill(HIST("mc_sec_pion_pt"), pt); + } + } + + if (passTpcKa && (pt < ptThreshold || passTofKa)) { + registryData.fill(HIST("rec_kaon_all"), pt); + registryData.fill(HIST("contamination_matrix_kaon"), pdg, pt); + if (isPrimary) { + if (pdg == 321) + registryData.fill(HIST("mc_rec_kaon_pt"), pt); + } else { + registryData.fill(HIST("mc_sec_kaon_pt"), pt); + } + } + + if (passTpcPr && (pt < ptThreshold || passTofPr)) { + registryData.fill(HIST("rec_proton_all"), pt); + registryData.fill(HIST("contamination_matrix_proton"), pdg, pt); + if (isPrimary) { + if (pdg == 2212) + registryData.fill(HIST("mc_rec_proton_pt"), pt); + } else { + registryData.fill(HIST("mc_sec_proton_pt"), pt); + } + } + } + } + PROCESS_SWITCH(jetHadronsPid, processMC, "Run on Monte Carlo", false); + + void processMCTruth(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles) + { + if (std::abs(mcCollision.posZ()) > zVtx) + return; + + for (auto const& mcpart : mcParticles) { + + if (!mcpart.isPhysicalPrimary()) + continue; + + if (mcpart.eta() < minEta || mcpart.eta() > maxEta) + continue; + if (mcpart.pt() < minPt || mcpart.pt() > maxPt) + continue; + + int pdg = std::abs(mcpart.pdgCode()); + double pt = mcpart.pt(); + + if (pdg == 211) { + registryData.fill(HIST("mc_gen_pion_pt"), pt); + } else if (pdg == 321) { + registryData.fill(HIST("mc_gen_kaon_pt"), pt); + } else if (pdg == 2212) { + registryData.fill(HIST("mc_gen_proton_pt"), pt); + } + } + } + PROCESS_SWITCH(jetHadronsPid, processMCTruth, "Run on Monte Carlo (Pure Truth)", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}