From 24dd1b212a06f45316aeaf7518bfd3c5888caef5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lars=20J=C3=B6rgensen?= Date: Thu, 27 Jun 2024 19:03:47 +0200 Subject: [PATCH 1/8] Add Run2 support for AngularCorrelationsInJets.cxx --- .../Nuspex/AngularCorrelationsInJets.cxx | 597 ++++++++++-------- 1 file changed, 343 insertions(+), 254 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx index ac3a77a04c2..ba24aa43e8e 100644 --- a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx +++ b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx @@ -8,24 +8,20 @@ // 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. -// -// Authors: Lars Jörgensen -// Date: 29.05.2024 -#include #include -#include -#include -#include "ReconstructionDataFormats/Track.h" +#include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "Framework/AnalysisDataModel.h" #include "Framework/ASoAHelpers.h" +#include "ReconstructionDataFormats/Track.h" +#include "CCDB/BasicCCDBManager.h" #include "Common/DataModel/TrackSelectionTables.h" #include "Common/DataModel/EventSelection.h" -#include "Framework/runDataProcessing.h" -#include "Framework/HistogramRegistry.h" -#include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/Core/PID/PIDTOF.h" +#include "Common/TableProducer/PID/pidTOFBase.h" #include "PWGLF/DataModel/LFParticleIdentification.h" #include "fastjet/PseudoJet.hh" #include "fastjet/AreaDefinition.hh" @@ -34,79 +30,21 @@ #include "fastjet/Selector.hh" #include "fastjet/tools/Subtractor.hh" #include "fastjet/tools/JetMedianBackgroundEstimator.hh" +// #include "PWGJE/Core/JetFinder.h" +// #include "PWGJE/Core/JetFindingUtilities.h" #include "TVector2.h" using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -// using SelectedCollisions = soa::Join; -using FullTracks = soa::Join; - struct AngularCorrelationsInJets { - - HistogramRegistry registryData{"jetOutput", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; - - void init(InitContext&) - { - // AxisSpec specName = {binningInfo, "axisLabel"}; - AxisSpec ptAxis = {1000, 0, 100, "#it{p}_{T} [GeV/#it{c}]"}; - AxisSpec particleTypeAxis = {4, 1, 5, "[p, ap, d, ad]"}; - AxisSpec nsigmapTAxis = {1000, -50, 50, "#it{p}_{T} [GeV/#it{c}]"}; - AxisSpec nsigmaAxis = {1000, -5, 5, "n#sigma"}; - AxisSpec dcazAxis = {200, -3, 3, "DCA_{z} [cm]"}; - AxisSpec dcaxyAxis = {200, -2, 2, "DCA_{xy} [cm]"}; - AxisSpec angDistPhiAxis = {1000, -2, 5, "#Delta#varphi"}; - AxisSpec angDistEtaAxis = {1000, -2, 2, "#Delta#eta"}; - - // Counters - registryData.add("hNumberOfEvents", "Number of events", HistType::kTH1I, {{1, 0, 1}}); - registryData.add("hNumberOfJets", "Total number of jets", HistType::kTH1I, {{1, 0, 1}}); - registryData.add("hEventProtocol", "Event protocol", HistType::kTH1I, {{20, 0, 20}}); - registryData.add("hTrackProtocol", "Track protocol", HistType::kTH1I, {{20, 0, 20}}); - registryData.add("hNumPartInJet", "Number of particles in a jet", HistType::kTH1I, {{200, 0, 200}}); - - // (Pseudo)Rapidity - registryData.add("hEtaFullEvent", "Particle pseudorapidity;#eta", HistType::kTH1F, {{200, -1, 1}}); - registryData.add("hJetRapidity", "Jet rapidity;#it{y}", HistType::kTH1F, {{200, -1, 1}}); - - // pT - registryData.add("hPtFullEvent", "p_{T} after basic cuts", HistType::kTH1F, {ptAxis}); - registryData.add("hPtJetParticle", "p_{T} of particles in jets", HistType::kTH1D, {ptAxis}); - registryData.add("hPtSubtractedJet", "Subtracted jet p_{T}", HistType::kTH1D, {ptAxis}); - registryData.add("hPtJetProtonDeuteron", "p_{T} of (anti)p, (anti)d", HistType::kTH2D, {particleTypeAxis, ptAxis}); - registryData.add("hPtTotalJet", "p_{T} of entire jet;#it{p}_{T} [GeV/#it{c}]", HistType::kTH1F, {{2000, 0, 500}}); - registryData.add("hPtDiff", "pT difference PseudoJet/original track;#it{p}_{T} [GeV/#it{c}]", HistType::kTH1D, {{100, -0.0000005, 0.0000005}}); - - // nSigma - registryData.add("hTPCnsigma", "TPC n#sigma for full event", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); - registryData.add("hTOFnsigma", "TOF n#sigma for full event", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); - registryData.add("hTPCnsigmaProton", "TPC n#sigma for (anti)proton", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); - registryData.add("hTOFnsigmaProton", "TOF n#sigma for (anti)proton", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); - registryData.add("hTPCnsigmaDeuteron", "TPC n#sigma for (anti)deuteron", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); - registryData.add("hTOFnsigmaDeuteron", "TOF n#sigma for (anti)deuteron", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); - - // DCA - registryData.add("hDCAxyFullJet", "DCA_{xy} of full jet", HistType::kTH2F, {ptAxis, dcaxyAxis}); - registryData.add("hDCAzFullJet", "DCA_{z} of full jet", HistType::kTH2F, {ptAxis, dcazAxis}); - registryData.add("hDCAzJetProton", "DCA_{z} of protons after TPC cut", HistType::kTH2F, {ptAxis, dcazAxis}); - registryData.add("hDCAzJetAntiproton", "DCA_{z} of antiprotons after TPC cut", HistType::kTH2F, {ptAxis, dcazAxis}); - registryData.add("hDCAzJetDeuteron", "DCA_{z} of deuterons after TPC cut", HistType::kTH2F, {ptAxis, dcazAxis}); - registryData.add("hDCAzJetAntideuteron", "DCA_{z} of antideuterons after TPC cut", HistType::kTH2F, {ptAxis, dcazAxis}); - - // Angular Distributions - registryData.add("hDeltaPhiSE", "#Delta#varphi of (anti)p, (anti)d in single event", HistType::kTH2D, {particleTypeAxis, angDistPhiAxis}); - registryData.add("hDeltaPhiME", "#Delta#varphi of (anti)p, (anti)d in mixed events", HistType::kTH2D, {particleTypeAxis, angDistPhiAxis}); - registryData.add("hDeltaPhiEtaSE", "#Delta#varphi vs #Delta#eta of (anti)p, (anti)d in single event", HistType::kTH3D, {particleTypeAxis, angDistPhiAxis, angDistEtaAxis}); - registryData.add("hDeltaPhiEtaME", "#Delta#varphi vs #Delta#eta of (anti)p, (anti)d in mixed events", HistType::kTH3D, {particleTypeAxis, angDistPhiAxis, angDistEtaAxis}); - - registryData.add("hJetConeRadius", "Jet Radius;#it{R}", HistType::kTH1F, {{100, 0, 1}}); - } - // Preliminary Cuts - // Configurable fTPCRefit{"TPCRefit", 0, "Require TPC refit"}; - Configurable fMinNCrossedRowsTPC{"minNCrossedRowsTPC", 70.0f, "min number of crossed rows TPC"}; - Configurable fMinReqClusterITS{"minReqClusterITS", 2.0, "min number of clusters required in ITS"}; + Configurable eventSelections{"eventSelection", "sel8", "choose event selection"}; + Configurable trackSelections{"trackSelection", "globalTracks", "set track selection"}; + Configurable fMinNCrossedRowsTPC{"minNCrossedRowsTPC", 70, "min number of crossed rows TPC"}; + Configurable fMinReqClusterITS{"minReqClusterITS", 2, "min number of clusters required in ITS"}; + Configurable fMinReqClusterTPC{"minReqClusterTPC", 70, "min number of clusters required in TPC"}; Configurable fMinRatioCrossedRowsTPC{"minRatioCrossedRowsTPC", 0.7f, "min ratio of crossed rows over findable clusters TPC"}; Configurable fMaxChi2ITS{"maxChi2ITS", 36.0f, "max chi2 per cluster ITS"}; Configurable fMaxChi2TPC{"maxChi2TPC", 4.0f, "max chi2 per cluster TPC"}; @@ -115,7 +53,7 @@ struct AngularCorrelationsInJets { Configurable fMaxEta{"maxEta", 0.8, "max pseudorapidity"}; // Jet Cuts - Configurable fJetRadius{"jetRadius", 0.4, "jet radius R"}; + Configurable fJetR{"jetR", 0.4, "jet resolution parameter"}; Configurable fMinJetPt{"minJetPt", 10.0, "minimum total pT to accept jet"}; Configurable fMinJetParticlePt{"minJetParticlePt", 0.0, "minimum pT to accept jet particle"}; Configurable fMinLeadingPt{"minLeadingPt", 5.0, "minimum pT for leading track"}; @@ -158,46 +96,133 @@ struct AngularCorrelationsInJets { Configurable fTrackBufferSize{"trackBufferSize", 2000, "Number of mixed-event tracks being stored"}; - std::vector fTrackBufferProton; - std::vector fTrackBufferAntiproton; - std::vector fTrackBufferDeuteron; - std::vector fTrackBufferAntideuteron; + Service ccdb; + // JetFinder jetFinder; + // std::vector constituents; + int mRunNumber; + + using FullTracksRun2 = soa::Join; + using FullTracksRun3 = soa::Join; + using BCsWithRun2Info = soa::Join; + + Filter prelimTrackCuts = (aod::track::itsChi2NCl < fMaxChi2ITS && + aod::track::tpcChi2NCl < fMaxChi2TPC && + nabs(aod::track::dcaXY) < fMaxDCAxy && + nabs(aod::track::dcaZ) < fMaxDCAz && + nabs(aod::track::eta) < fMaxEta); - //**************************************************************************************************** + Preslice perCollisionFullTracksRun2 = o2::aod::track::collisionId; - template - void angCorrData(const CollisionType& /*collision*/, const TracksType& tracks) + AxisSpec ptAxis = {1000, 0, 100, "#it{p}_{T} [GeV/#it{c}]"}; + AxisSpec particleTypeAxis = {4, 1, 5, "[p, ap, d, ad]"}; + AxisSpec nsigmapTAxis = {1000, -50, 50, "#it{p}_{T} [GeV/#it{c}]"}; + AxisSpec nsigmaAxis = {1000, -5, 5, "n#sigma"}; + AxisSpec dcazAxis = {200, -3, 3, "DCA_{z} [cm]"}; + AxisSpec dcaxyAxis = {200, -2, 2, "DCA_{xy} [cm]"}; + AxisSpec angDistPhiAxis = {1000, -2, 5, "#Delta#varphi"}; + AxisSpec angDistEtaAxis = {1000, -2, 2, "#Delta#eta"}; + + HistogramRegistry registryData{"dataOutput", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + HistogramRegistry registryQA{"dataQA", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; + + void init(o2::framework::InitContext&) { - registryData.fill(HIST("hNumberOfEvents"), 0); - registryData.fill(HIST("hEventProtocol"), 0); + mRunNumber = 0; + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); - std::vector jetInput; - std::vector jets; - std::vector constituents; - fastjet::PseudoJet hardestJet(0., 0., 0., 0.); - fastjet::PseudoJet subtractedJet(0., 0., 0., 0.); - jetInput.clear(); - jets.clear(); - constituents.clear(); + // Counters + registryData.add("hNumberOfEvents", "Number of events", HistType::kTH1I, {{1, 0, 1}}); + registryData.add("hNumberOfJets", "Total number of jets", HistType::kTH1I, {{1, 0, 1}}); + registryData.add("hEventProtocol", "Event protocol", HistType::kTH1I, {{20, 0, 20}}); + registryData.add("hTrackProtocol", "Track protocol", HistType::kTH1I, {{20, 0, 20}}); + registryData.add("hNumPartInJet", "Number of particles in a jet", HistType::kTH1I, {{200, 0, 200}}); - for (const auto& track : tracks) { - registryData.fill(HIST("hTrackProtocol"), 0); - registryData.fill(HIST("hEtaFullEvent"), track.eta()); - registryData.fill(HIST("hPtFullEvent"), track.pt()); + // (Pseudo)Rapidity + registryData.add("hEtaFullEvent", "Particle pseudorapidity;#eta", HistType::kTH1F, {{200, -1, 1}}); + registryData.add("hJetRapidity", "Jet rapidity;#it{y}", HistType::kTH1F, {{200, -1, 1}}); - fastjet::PseudoJet inputPseudoJet(track.px(), track.py(), track.pz(), track.energy(track.mass())); // is this fine? just TOF for invariant mass? - inputPseudoJet.set_user_index(track.globalIndex()); - jetInput.emplace_back(inputPseudoJet); - } + // pT + registryData.add("hPtFullEvent", "p_{T} after basic cuts", HistType::kTH1F, {ptAxis}); + registryData.add("hPtJetParticle", "p_{T} of particles in jets", HistType::kTH1D, {ptAxis}); + registryData.add("hPtSubtractedJet", "Subtracted jet p_{T}", HistType::kTH1D, {ptAxis}); + registryData.add("hPtJetProtonDeuteron", "p_{T} of (anti)p, (anti)d", HistType::kTH2D, {particleTypeAxis, ptAxis}); + registryData.add("hPtTotalJet", "p_{T} of entire jet;#it{p}_{T} [GeV/#it{c}]", HistType::kTH1F, {{2000, 0, 500}}); + registryData.add("hPtDiff", "pT difference PseudoJet/original track;#it{p}_{T} [GeV/#it{c}]", HistType::kTH1D, {{100, -5, 5}}); - if (jetInput.size() < 2) + // nSigma + registryData.add("hTPCsignal", "TPC signal", HistType::kTH2F, {nsigmapTAxis, {1000, -100, 100}}); + registryData.add("hTOFsignal", "TOF signal", HistType::kTH2F, {nsigmapTAxis, {200, -1, 1}}); + registryData.add("hTPCnsigmaProton", "TPC n#sigma for (anti)proton", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); + registryData.add("hTOFnsigmaProton", "TOF n#sigma for (anti)proton", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); + registryData.add("hTPCnsigmaDeuteron", "TPC n#sigma for (anti)deuteron", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); + registryData.add("hTOFnsigmaDeuteron", "TOF n#sigma for (anti)deuteron", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); + + // DCA + registryData.add("hDCAxyFullJet", "DCA_{xy} of full jet", HistType::kTH2F, {ptAxis, dcaxyAxis}); + registryData.add("hDCAzFullJet", "DCA_{z} of full jet", HistType::kTH2F, {ptAxis, dcazAxis}); + registryData.add("hDCAzJetProton", "DCA_{z} of protons after TPC cut", HistType::kTH2F, {ptAxis, dcazAxis}); + registryData.add("hDCAzJetAntiproton", "DCA_{z} of antiprotons after TPC cut", HistType::kTH2F, {ptAxis, dcazAxis}); + registryData.add("hDCAzJetDeuteron", "DCA_{z} of deuterons after TPC cut", HistType::kTH2F, {ptAxis, dcazAxis}); + registryData.add("hDCAzJetAntideuteron", "DCA_{z} of antideuterons after TPC cut", HistType::kTH2F, {ptAxis, dcazAxis}); + + // Angular Distributions + registryData.add("hDeltaPhiSE", "#Delta#varphi of (anti)p, (anti)d in single event", HistType::kTH2D, {particleTypeAxis, angDistPhiAxis}); + registryData.add("hDeltaPhiME", "#Delta#varphi of (anti)p, (anti)d in mixed events", HistType::kTH2D, {particleTypeAxis, angDistPhiAxis}); + registryData.add("hDeltaPhiEtaSE", "#Delta#varphi vs #Delta#eta of (anti)p, (anti)d in single event", HistType::kTH3D, {particleTypeAxis, angDistPhiAxis, angDistEtaAxis}); + registryData.add("hDeltaPhiEtaME", "#Delta#varphi vs #Delta#eta of (anti)p, (anti)d in mixed events", HistType::kTH3D, {particleTypeAxis, angDistPhiAxis, angDistEtaAxis}); + + registryData.add("hJetConeRadius", "Jet Radius;#it{R}", HistType::kTH1F, {{100, 0, 1}}); + + //TODO: add QA histograms + } + + std::vector fTrackBufferProton; + std::vector fTrackBufferAntiproton; + std::vector fTrackBufferDeuteron; + std::vector fTrackBufferAntideuteron; + //TODO: check if FullTracksRun2 works for Run3 too or add Run3 track buffers + + template + void initCCDB(Bc const& bc) + { + if (mRunNumber == bc.runNumber()) { return; - registryData.fill(HIST("hEventProtocol"), 1); + } + mRunNumber = bc.runNumber(); + } + + template + bool selectTrack(T const& track) + { + if (track.tpcNClsCrossedRows() < fMinRatioCrossedRowsTPC * track.tpcNClsFindable() || + track.tpcNClsCrossedRows() < fMinNCrossedRowsTPC || + track.tpcNClsFound() < fMinReqClusterTPC || + track.itsNCls() < fMinReqClusterITS) { + return false; + } + if (doprocessRun2) { + if (!(track.trackType() & o2::aod::track::Run2Track)/* || + !(track.flags() & o2::aod::track::TPCrefit) || + !(track.flags() & o2::aod::track::ITSrefit) */) { + return false; + } + } + return true; + } + + void fillConstituents(std::vector jetInput, std::vector constituents) { + std::vector jets; + jets.clear(); + fastjet::PseudoJet hardestJet(0., 0., 0., 0.); + fastjet::PseudoJet subtractedJet(0., 0., 0., 0.); double ghost_maxrap = 1.0; double ghost_area = 0.005; int ghost_repeat = 1; - fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, fJetRadius); + fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, fJetR); fastjet::JetDefinition jetDefBkg(fastjet::kt_algorithm, 0.5); fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(ghost_maxrap, ghost_repeat, ghost_area)); fastjet::AreaDefinition areaDefBkg(fastjet::active_area_explicit_ghosts, fastjet::GhostedAreaSpec(ghost_maxrap)); @@ -205,27 +230,27 @@ struct AngularCorrelationsInJets { jets = sorted_by_pt(clusterSeq.inclusive_jets()); if (jets.size() == 0) return; - registryData.fill(HIST("hEventProtocol"), 2); + registryData.fill(HIST("hEventProtocol"), 3); hardestJet = jets[0]; if (hardestJet.pt() < fMinJetPt) return; - registryData.fill(HIST("hEventProtocol"), 3); + registryData.fill(HIST("hEventProtocol"), 4); if (hardestJet.constituents().size() < 2) return; - registryData.fill(HIST("hEventProtocol"), 4); + registryData.fill(HIST("hEventProtocol"), 5); registryData.fill(HIST("hNumberOfJets"), 0); registryData.fill(HIST("hPtTotalJet"), hardestJet.pt()); registryData.fill(HIST("hJetRapidity"), hardestJet.rap()); constituents = hardestJet.constituents(); - for (int i = 0; i < constituents.size(); i++) { + for (int i = 0; i < static_cast(constituents.size()); i++) { registryData.fill(HIST("hPtJetParticle"), constituents[i].pt()); double DeltaPhi = TVector2::Phi_0_2pi(constituents[i].phi() - hardestJet.phi()); double DeltaEta = constituents[i].eta() - hardestJet.eta(); - double Delta = TMath::Sqrt(DeltaPhi * DeltaPhi + DeltaEta * DeltaEta) / (constituents[i].pt() * constituents[i].pt()); // need 1/pT^2? + double Delta = TMath::Sqrt(DeltaPhi * DeltaPhi + DeltaEta * DeltaEta); // need 1/pT^2? registryData.fill(HIST("hJetConeRadius"), Delta); } @@ -237,127 +262,14 @@ struct AngularCorrelationsInJets { subtractedJet = subtractor(hardestJet); if (subtractedJet.has_constituents()) { - for (int i = 0; i < subtractedJet.constituents().size(); i++) { + for (int i = 0; i < static_cast(subtractedJet.constituents().size()); i++) { registryData.fill(HIST("hPtSubtractedJet"), subtractedJet.constituents()[i].pt()); } } - - std::vector jetProtons; - std::vector jetAntiprotons; - std::vector jetDeuterons; - std::vector jetAntideuterons; - - for (int i = 0; i < constituents.size(); i++) { - fastjet::PseudoJet pseudoParticle = constituents[i]; - int id = pseudoParticle.user_index(); - typename TracksType::iterator jetParticle = tracks.iteratorAt(id); - registryData.fill(HIST("hDCAxyFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaXY()); - registryData.fill(HIST("hDCAzFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaZ()); - registryData.fill(HIST("hTPCnsigma"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcNSigmaPr()); - registryData.fill(HIST("hTOFnsigma"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaPr()); - - double ptDiff = pseudoParticle.pt() - jetParticle.pt(); - registryData.fill(HIST("hPtDiff"), ptDiff); - int particleType = 0; - - if (jetParticle.pt() < fMinJetParticlePt) - continue; - if (IsProton(jetParticle) || IsAntiproton(jetParticle)) { // collect (anti)protons in jet - registryData.fill(HIST("hTPCnsigmaProton"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcNSigmaPr()); - registryData.fill(HIST("hTOFnsigmaProton"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaPr()); - if (IsProton(jetParticle)) { - particleType = 1; - registryData.fill(HIST("hTrackProtocol"), 6); // # protons - jetProtons.emplace_back(jetParticle); - } else { - particleType = 2; - registryData.fill(HIST("hTrackProtocol"), 7); // # antiprotons - jetAntiprotons.emplace_back(jetParticle); - } - } else if (IsDeuteron(jetParticle) || IsAntideuteron(jetParticle)) { // collect (anti)deuterons in jet - registryData.fill(HIST("hTPCnsigmaDeuteron"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcNSigmaDe()); - registryData.fill(HIST("hTOFnsigmaDeuteron"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaDe()); - if (IsDeuteron(jetParticle)) { - particleType = 3; - registryData.fill(HIST("hTrackProtocol"), 8); // # deuterons - jetDeuterons.emplace_back(jetParticle); - } else { - particleType = 4; - registryData.fill(HIST("hTrackProtocol"), 9); // # antideuterons - jetAntideuterons.emplace_back(jetParticle); - } - } - registryData.fill(HIST("hPtJetProtonDeuteron"), particleType, jetParticle.pt()); - } // for (int i=0; i<(int)constituents.size(); i++) - - if ((jetProtons.size() < 2) && (jetAntiprotons.size() < 2) && (jetDeuterons.size() < 2) && (jetAntideuterons.size() < 2)) - return; - registryData.fill(HIST("hEventProtocol"), 5); - - if (jetProtons.size() > 1) { - for (int i = 0; i < jetProtons.size(); i++) { - for (int j = i + 1; j < jetProtons.size(); j++) { - double DeltaPhi = TVector2::Phi_0_2pi(jetProtons[i].phi() - jetProtons[j].phi()); - double DeltaEta = TMath::Abs(jetProtons[i].eta() - jetProtons[j].eta()); - if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); - } - registryData.fill(HIST("hDeltaPhiSE"), 1, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaSE"), 1, DeltaPhi, DeltaEta); - } - FillMixedEventDeltas(jetProtons[i], 1); - SetTrackBuffer(jetProtons[i], 1); - } - } - if (jetAntiprotons.size() > 1) { - for (int i = 0; i < jetAntiprotons.size(); i++) { - for (int j = i + 1; j < jetAntiprotons.size(); j++) { - double DeltaPhi = TVector2::Phi_0_2pi(jetAntiprotons[i].phi() - jetAntiprotons[j].phi()); - double DeltaEta = TMath::Abs(jetAntiprotons[i].eta() - jetAntiprotons[j].eta()); - if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); - } - registryData.fill(HIST("hDeltaPhiSE"), 2, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaSE"), 2, DeltaPhi, DeltaEta); - } - FillMixedEventDeltas(jetAntiprotons[i], 2); - SetTrackBuffer(jetAntiprotons[i], 2); - } - } - if (jetDeuterons.size() > 1) { - for (int i = 0; i < jetDeuterons.size(); i++) { - for (int j = i + 1; j < jetDeuterons.size(); j++) { - double DeltaPhi = TVector2::Phi_0_2pi(jetDeuterons[i].phi() - jetDeuterons[j].phi()); - double DeltaEta = TMath::Abs(jetDeuterons[i].eta() - jetDeuterons[j].eta()); - if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); - } - registryData.fill(HIST("hDeltaPhiSE"), 3, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaSE"), 3, DeltaPhi, DeltaEta); - } - FillMixedEventDeltas(jetDeuterons[i], 3); - SetTrackBuffer(jetDeuterons[i], 3); - } - } - if (jetAntideuterons.size() > 1) { - for (int i = 0; i < jetAntideuterons.size(); i++) { - for (int j = i + 1; j < jetAntideuterons.size(); j++) { - double DeltaPhi = TVector2::Phi_0_2pi(jetAntideuterons[i].phi() - jetAntideuterons[j].phi()); - double DeltaEta = TMath::Abs(jetAntideuterons[i].eta() - jetAntideuterons[j].eta()); - if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); - } - registryData.fill(HIST("hDeltaPhiSE"), 4, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaSE"), 4, DeltaPhi, DeltaEta); - } - FillMixedEventDeltas(jetAntideuterons[i], 4); - SetTrackBuffer(jetAntideuterons[i], 4); - } - } } - template - bool IsProton(const T1& track) + template + bool isProton(const T& track) { bool isProton = false; if (track.sign() < 0) @@ -383,8 +295,8 @@ struct AngularCorrelationsInJets { return isProton; } - template - bool IsAntiproton(const T2& track) + template + bool isAntiproton(const T& track) { bool isAntiproton = false; if (track.sign() < 0) @@ -410,8 +322,8 @@ struct AngularCorrelationsInJets { return isAntiproton; } - template - bool IsDeuteron(const T3& track) + template + bool isDeuteron(const T& track) { bool isDeuteron = false; if (track.sign() < 0) @@ -437,8 +349,8 @@ struct AngularCorrelationsInJets { return isDeuteron; } - template - bool IsAntideuteron(const T4& track) + template + bool isAntideuteron(const T& track) { bool isAntideuteron = false; if (track.sign() < 0) @@ -464,39 +376,39 @@ struct AngularCorrelationsInJets { return isAntideuteron; } - template - void SetTrackBuffer(const T5& track, int particleType) + template + void setTrackBuffer(const T& track, int particleType) //TODO: buffer as argument, lose switch statement { switch (particleType) { case 1: - if (fTrackBufferProton.size() == fTrackBufferSize) { + if (static_cast(fTrackBufferProton.size()) == fTrackBufferSize) { fTrackBufferProton.insert(fTrackBufferProton.begin(), track); fTrackBufferProton.resize(fTrackBufferSize); - } else if (fTrackBufferProton.size() < fTrackBufferSize) { + } else if (static_cast(fTrackBufferProton.size()) < fTrackBufferSize) { fTrackBufferProton.emplace_back(track); } break; case 2: - if (fTrackBufferAntiproton.size() == fTrackBufferSize) { + if (static_cast(fTrackBufferAntiproton.size()) == fTrackBufferSize) { fTrackBufferAntiproton.insert(fTrackBufferAntiproton.begin(), track); fTrackBufferAntiproton.resize(fTrackBufferSize); - } else if (fTrackBufferAntiproton.size() < fTrackBufferSize) { + } else if (static_cast(fTrackBufferAntiproton.size()) < fTrackBufferSize) { fTrackBufferAntiproton.emplace_back(track); } break; case 3: - if (fTrackBufferDeuteron.size() == fTrackBufferSize) { + if (static_cast(fTrackBufferDeuteron.size()) == fTrackBufferSize) { fTrackBufferDeuteron.insert(fTrackBufferDeuteron.begin(), track); fTrackBufferDeuteron.resize(fTrackBufferSize); - } else if (fTrackBufferDeuteron.size() < fTrackBufferSize) { + } else if (static_cast(fTrackBufferDeuteron.size()) < fTrackBufferSize) { fTrackBufferDeuteron.emplace_back(track); } break; case 4: - if (fTrackBufferAntideuteron.size() == fTrackBufferSize) { + if (static_cast(fTrackBufferAntideuteron.size()) == fTrackBufferSize) { fTrackBufferAntideuteron.insert(fTrackBufferAntideuteron.begin(), track); fTrackBufferAntideuteron.resize(fTrackBufferSize); - } else if (fTrackBufferAntideuteron.size() < fTrackBufferSize) { + } else if (static_cast(fTrackBufferAntideuteron.size()) < fTrackBufferSize) { fTrackBufferAntideuteron.emplace_back(track); } break; @@ -505,14 +417,14 @@ struct AngularCorrelationsInJets { } } - template - void FillMixedEventDeltas(const T6& track, int particleType) + template + void fillMixedEventDeltas(const T& track, int particleType) //TODO: buffer as argument, lose switch statement { switch (particleType) { case 1: if (fTrackBufferProton.size() == 0) return; - for (int i = 0; i < fTrackBufferProton.size(); i++) { // can I do this even if the track buffer isn't even full yet? + for (int i = 0; i < static_cast(fTrackBufferProton.size()); i++) { // even if the track buffer isn't even full yet? double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - fTrackBufferProton[i].phi()); double DeltaEta = TMath::Abs(track.eta() - fTrackBufferProton[i].eta()); registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); @@ -522,7 +434,7 @@ struct AngularCorrelationsInJets { case 2: if (fTrackBufferAntiproton.size() == 0) return; - for (int i = 0; i < fTrackBufferAntiproton.size(); i++) { + for (int i = 0; i < static_cast(fTrackBufferAntiproton.size()); i++) { double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - fTrackBufferAntiproton[i].phi()); double DeltaEta = TMath::Abs(track.eta() - fTrackBufferAntiproton[i].eta()); registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); @@ -532,7 +444,7 @@ struct AngularCorrelationsInJets { case 3: if (fTrackBufferDeuteron.size() == 0) return; - for (int i = 0; i < fTrackBufferDeuteron.size(); i++) { + for (int i = 0; i < static_cast(fTrackBufferDeuteron.size()); i++) { double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - fTrackBufferDeuteron[i].phi()); double DeltaEta = TMath::Abs(track.eta() - fTrackBufferDeuteron[i].eta()); registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); @@ -542,7 +454,7 @@ struct AngularCorrelationsInJets { case 4: if (fTrackBufferAntideuteron.size() == 0) return; - for (int i = 0; i < fTrackBufferAntideuteron.size(); i++) { + for (int i = 0; i < static_cast(fTrackBufferAntideuteron.size()); i++) { double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - fTrackBufferAntideuteron[i].phi()); double DeltaEta = TMath::Abs(track.eta() - fTrackBufferAntideuteron[i].eta()); registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); @@ -554,21 +466,198 @@ struct AngularCorrelationsInJets { } } - //**************************************************************************************************** + template + void fillHistogramsRun2(T const& collision, U const& allTracks) { + std::vector jetInput; + std::vector constituents; + jetInput.clear(); + constituents.clear(); + + auto tracks = allTracks.sliceBy(perCollisionFullTracksRun2, collision.globalIndex()); + + for (const auto& track : tracks) { + if (!selectTrack(track)) + continue; + + registryData.fill(HIST("hPtFullEvent"), track.pt()); + fastjet::PseudoJet inputPseudoJet(track.px(), track.py(), track.pz(), track.energy(track.mass())); + inputPseudoJet.set_user_index(track.globalIndex()+1); + jetInput.emplace_back(inputPseudoJet); + } + + if (jetInput.size() < 2) + return; + registryData.fill(HIST("hEventProtocol"), 2); + + fillConstituents(jetInput, constituents); + + std::vector jetProtons; + std::vector jetAntiprotons; + std::vector jetDeuterons; + std::vector jetAntideuterons; + + for (int i = 0; i < static_cast(constituents.size()); i++) { + fastjet::PseudoJet pseudoParticle = constituents[i]; + uint64_t id = pseudoParticle.user_index(); + auto jetParticle = tracks.iteratorAt(id); + registryData.fill(HIST("hDCAxyFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaXY()); + registryData.fill(HIST("hDCAzFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaZ()); + registryData.fill(HIST("hTPCsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcSignal()); + registryData.fill(HIST("hTOFsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofSignal()); + + double ptDiff = pseudoParticle.pt() - jetParticle.pt(); + registryData.fill(HIST("hPtDiff"), ptDiff); + int particleType = 0; + + if (jetParticle.pt() < fMinJetParticlePt) + continue; + if (isProton(jetParticle) || isAntiproton(jetParticle)) { // collect (anti)protons in jet + registryData.fill(HIST("hTPCnsigmaProton"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcNSigmaPr()); + registryData.fill(HIST("hTOFnsigmaProton"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaPr()); + if (isProton(jetParticle)) { + particleType = 1; + registryData.fill(HIST("hTrackProtocol"), 6); // # protons + jetProtons.emplace_back(jetParticle); + } else { + particleType = 2; + registryData.fill(HIST("hTrackProtocol"), 7); // # antiprotons + jetAntiprotons.emplace_back(jetParticle); + } + } else if (isDeuteron(jetParticle) || isAntideuteron(jetParticle)) { // collect (anti)deuterons in jet + registryData.fill(HIST("hTPCnsigmaDeuteron"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcNSigmaDe()); + registryData.fill(HIST("hTOFnsigmaDeuteron"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaDe()); + if (isDeuteron(jetParticle)) { + particleType = 3; + registryData.fill(HIST("hTrackProtocol"), 8); // # deuterons + jetDeuterons.emplace_back(jetParticle); + } else { + particleType = 4; + registryData.fill(HIST("hTrackProtocol"), 9); // # antideuterons + jetAntideuterons.emplace_back(jetParticle); + } + } + registryData.fill(HIST("hPtJetProtonDeuteron"), particleType, jetParticle.pt()); + } // for (int i=0; i(constituents.size()); i++) + + if ((jetProtons.size() < 2) && (jetAntiprotons.size() < 2) && (jetDeuterons.size() < 2) && (jetAntideuterons.size() < 2)) + return; + registryData.fill(HIST("hEventProtocol"), 5); + + if (jetProtons.size() > 1) { //TODO: write function that takes jetParticle vector, lose 3/4 if statements + for (int i = 0; i < static_cast(jetProtons.size()); i++) { + for (int j = i + 1; j < static_cast(jetProtons.size()); j++) { + double DeltaPhi = TVector2::Phi_0_2pi(jetProtons[i].phi() - jetProtons[j].phi()); + double DeltaEta = TMath::Abs(jetProtons[i].eta() - jetProtons[j].eta()); + if (DeltaPhi > (1.5 * TMath::Pi())) { + DeltaPhi = DeltaPhi - 2 * TMath::Pi(); + } + registryData.fill(HIST("hDeltaPhiSE"), 1, DeltaPhi); + registryData.fill(HIST("hDeltaPhiEtaSE"), 1, DeltaPhi, DeltaEta); + } + fillMixedEventDeltas(jetProtons[i], 1); + setTrackBuffer(jetProtons[i], 1); + } + } + if (jetAntiprotons.size() > 1) { + for (int i = 0; i < static_cast(jetAntiprotons.size()); i++) { + for (int j = i + 1; j < static_cast(jetAntiprotons.size()); j++) { + double DeltaPhi = TVector2::Phi_0_2pi(jetAntiprotons[i].phi() - jetAntiprotons[j].phi()); + double DeltaEta = TMath::Abs(jetAntiprotons[i].eta() - jetAntiprotons[j].eta()); + if (DeltaPhi > (1.5 * TMath::Pi())) { + DeltaPhi = DeltaPhi - 2 * TMath::Pi(); + } + registryData.fill(HIST("hDeltaPhiSE"), 2, DeltaPhi); + registryData.fill(HIST("hDeltaPhiEtaSE"), 2, DeltaPhi, DeltaEta); + } + fillMixedEventDeltas(jetAntiprotons[i], 2); + setTrackBuffer(jetAntiprotons[i], 2); + } + } + if (jetDeuterons.size() > 1) { + for (int i = 0; i < static_cast(jetDeuterons.size()); i++) { + for (int j = i + 1; j < static_cast(jetDeuterons.size()); j++) { + double DeltaPhi = TVector2::Phi_0_2pi(jetDeuterons[i].phi() - jetDeuterons[j].phi()); + double DeltaEta = TMath::Abs(jetDeuterons[i].eta() - jetDeuterons[j].eta()); + if (DeltaPhi > (1.5 * TMath::Pi())) { + DeltaPhi = DeltaPhi - 2 * TMath::Pi(); + } + registryData.fill(HIST("hDeltaPhiSE"), 3, DeltaPhi); + registryData.fill(HIST("hDeltaPhiEtaSE"), 3, DeltaPhi, DeltaEta); + } + fillMixedEventDeltas(jetDeuterons[i], 3); + setTrackBuffer(jetDeuterons[i], 3); + } + } + if (jetAntideuterons.size() > 1) { + for (int i = 0; i < static_cast(jetAntideuterons.size()); i++) { + for (int j = i + 1; j < static_cast(jetAntideuterons.size()); j++) { + double DeltaPhi = TVector2::Phi_0_2pi(jetAntideuterons[i].phi() - jetAntideuterons[j].phi()); + double DeltaEta = TMath::Abs(jetAntideuterons[i].eta() - jetAntideuterons[j].eta()); + if (DeltaPhi > (1.5 * TMath::Pi())) { + DeltaPhi = DeltaPhi - 2 * TMath::Pi(); + } + registryData.fill(HIST("hDeltaPhiSE"), 4, DeltaPhi); + registryData.fill(HIST("hDeltaPhiEtaSE"), 4, DeltaPhi, DeltaEta); + } + fillMixedEventDeltas(jetAntideuterons[i], 4); + setTrackBuffer(jetAntideuterons[i], 4); + } + } + } + + template + void fillHistogramsRun3(T const& collision, U const& tracks) { + std::vector jetInput; + std::vector constituents; + jetInput.clear(); + constituents.clear(); - Filter prelimTrackCuts = (/* aod::track::TPCrefit == fTPCRefit && */ aod::track::itsChi2NCl < fMaxChi2ITS && aod::track::tpcChi2NCl < fMaxChi2TPC && nabs(aod::track::dcaXY) < fMaxDCAxy && nabs(aod::track::dcaZ) < fMaxDCAz && nabs(aod::track::eta) < fMaxEta); + for (const auto& track : tracks) { + registryData.fill(HIST("hTrackProtocol"), 0); + registryData.fill(HIST("hEtaFullEvent"), track.eta()); + registryData.fill(HIST("hPtFullEvent"), track.pt()); - void process_ang_corr_data(aod::Collision const& collision, soa::Filtered const& tracks) + fastjet::PseudoJet inputPseudoJet(track.px(), track.py(), track.pz(), track.energy(track.mass())); + inputPseudoJet.set_user_index(track.globalIndex()); + jetInput.emplace_back(inputPseudoJet); + } + + fillConstituents(jetInput, constituents); + //TODO: add rest of the logic when it is are optimised + } + + void processRun2(soa::Join const& collisions, + soa::Filtered const& tracks, + BCsWithRun2Info const&) + { + int index = 0; + for (const auto& collision : collisions) { + index++; + auto bc = collision.bc_as(); + initCCDB(bc); + + registryData.fill(HIST("hEventProtocol"), 0); + registryData.fill(HIST("hNumberOfEvents"), 0); + if (!collision.alias_bit(kINT7)) + continue; + registryData.fill(HIST("hEventProtocol"), 1); + + fillHistogramsRun2(collision, tracks); + } + } + PROCESS_SWITCH(AngularCorrelationsInJets, processRun2, "process Run 2 data", true); + + void processRun3(aod::Collision const& collision, soa::Filtered const& tracks) { - angCorrData(collision, tracks); + fillHistogramsRun3(collision, tracks); } - PROCESS_SWITCH(AngularCorrelationsInJets, process_ang_corr_data, "ang correlations in reco'ed jets", true); + PROCESS_SWITCH(AngularCorrelationsInJets, processRun3, "process Run 3 data", false); }; -//**************************************************************************************************** + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc, TaskName{"angular-correlations-in-jets"})}; + adaptAnalysisTask(cfgc)}; } From dd6e16af65a20b07072cb530515b62ffc481824f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lars=20J=C3=B6rgensen?= Date: Sun, 30 Jun 2024 17:19:31 +0200 Subject: [PATCH 2/8] AngularCorrelationsInJets.cxx major overhaul, reduced line count --- .../Nuspex/AngularCorrelationsInJets.cxx | 284 ++++++++---------- 1 file changed, 120 insertions(+), 164 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx index ba24aa43e8e..2ccf57883eb 100644 --- a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx +++ b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx @@ -8,6 +8,8 @@ // 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. +// +/// \author Lars Jörgensen #include @@ -30,9 +32,9 @@ #include "fastjet/Selector.hh" #include "fastjet/tools/Subtractor.hh" #include "fastjet/tools/JetMedianBackgroundEstimator.hh" +#include "TVector2.h" // #include "PWGJE/Core/JetFinder.h" // #include "PWGJE/Core/JetFindingUtilities.h" -#include "TVector2.h" using namespace o2; using namespace o2::framework; @@ -153,8 +155,8 @@ struct AngularCorrelationsInJets { registryData.add("hPtDiff", "pT difference PseudoJet/original track;#it{p}_{T} [GeV/#it{c}]", HistType::kTH1D, {{100, -5, 5}}); // nSigma - registryData.add("hTPCsignal", "TPC signal", HistType::kTH2F, {nsigmapTAxis, {1000, -100, 100}}); - registryData.add("hTOFsignal", "TOF signal", HistType::kTH2F, {nsigmapTAxis, {200, -1, 1}}); + registryData.add("hTPCsignal", "TPC signal", HistType::kTH2F, {{1000, 0, 100, "#it{p} [GeV/#it{c}]"}, {5000, 0, 5000, "d#it{E}/d#it{X} (a.u.)"}}); + registryData.add("hTOFsignal", "TOF signal", HistType::kTH2F, {{1000, 0, 100, "#it{p} [GeV/#it{c}]"}, {550, 0, 1.1, "#beta (TOF)"}}); registryData.add("hTPCnsigmaProton", "TPC n#sigma for (anti)proton", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); registryData.add("hTOFnsigmaProton", "TOF n#sigma for (anti)proton", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); registryData.add("hTPCnsigmaDeuteron", "TPC n#sigma for (anti)deuteron", HistType::kTH2F, {nsigmapTAxis, nsigmaAxis}); @@ -177,6 +179,9 @@ struct AngularCorrelationsInJets { registryData.add("hJetConeRadius", "Jet Radius;#it{R}", HistType::kTH1F, {{100, 0, 1}}); //TODO: add QA histograms + registryQA.add("hTOFmass", "TOF mass", HistType::kTH2F, {ptAxis, {1000, 0, 5, "#it{m} [GeV/#it{c}^{2}]"}}); + registryQA.get(HIST("hTOFmass"))->Sumw2(); + } std::vector fTrackBufferProton; @@ -213,9 +218,12 @@ struct AngularCorrelationsInJets { return true; } - void fillConstituents(std::vector jetInput, std::vector constituents) { + std::vector findJets(std::vector jetInput) + { std::vector jets; + std::vector constituents; jets.clear(); + constituents.clear(); fastjet::PseudoJet hardestJet(0., 0., 0., 0.); fastjet::PseudoJet subtractedJet(0., 0., 0., 0.); @@ -229,27 +237,26 @@ struct AngularCorrelationsInJets { fastjet::ClusterSequenceArea clusterSeq(jetInput, jetDef, areaDef); // or CSActiveArea? jets = sorted_by_pt(clusterSeq.inclusive_jets()); if (jets.size() == 0) - return; + return constituents; + registryData.fill(HIST("hEventProtocol"), 3); hardestJet = jets[0]; if (hardestJet.pt() < fMinJetPt) - return; + return constituents; registryData.fill(HIST("hEventProtocol"), 4); if (hardestJet.constituents().size() < 2) - return; + return constituents; registryData.fill(HIST("hEventProtocol"), 5); registryData.fill(HIST("hNumberOfJets"), 0); registryData.fill(HIST("hPtTotalJet"), hardestJet.pt()); registryData.fill(HIST("hJetRapidity"), hardestJet.rap()); - constituents = hardestJet.constituents(); - - for (int i = 0; i < static_cast(constituents.size()); i++) { - registryData.fill(HIST("hPtJetParticle"), constituents[i].pt()); - double DeltaPhi = TVector2::Phi_0_2pi(constituents[i].phi() - hardestJet.phi()); - double DeltaEta = constituents[i].eta() - hardestJet.eta(); + for (const auto& constituent : hardestJet.constituents()) { + registryData.fill(HIST("hPtJetParticle"), constituent.pt()); + double DeltaPhi = TVector2::Phi_0_2pi(constituent.phi() - hardestJet.phi()); + double DeltaEta = constituent.eta() - hardestJet.eta(); double Delta = TMath::Sqrt(DeltaPhi * DeltaPhi + DeltaEta * DeltaEta); // need 1/pT^2? registryData.fill(HIST("hJetConeRadius"), Delta); } @@ -262,10 +269,11 @@ struct AngularCorrelationsInJets { subtractedJet = subtractor(hardestJet); if (subtractedJet.has_constituents()) { - for (int i = 0; i < static_cast(subtractedJet.constituents().size()); i++) { - registryData.fill(HIST("hPtSubtractedJet"), subtractedJet.constituents()[i].pt()); + for (const auto& subConstituent : subtractedJet.constituents()) { + registryData.fill(HIST("hPtSubtractedJet"), subConstituent.pt()); } } + return hardestJet.constituents(); } template @@ -275,6 +283,7 @@ struct AngularCorrelationsInJets { if (track.sign() < 0) return isProton; + //TPC if (track.pt() < fProtonTPCTOFpT && TMath::Abs(track.tpcNSigmaPr()) < fProtonTPCnsigLow) isProton = true; if (track.pt() > fProtonTPCTOFpT && TMath::Abs(track.tpcNSigmaPr()) < fProtonTPCnsigHigh) @@ -287,6 +296,7 @@ struct AngularCorrelationsInJets { if (TMath::Abs(track.dcaZ()) > fProtonDCAz) return false; + //TOF if (track.pt() < fProtonTPCTOFpT && TMath::Abs(track.tofNSigmaPr()) < fProtonTOFnsigLow) isProton = true; if (track.pt() > fProtonTPCTOFpT && TMath::Abs(track.tofNSigmaPr()) < fProtonTOFnsigHigh) @@ -302,6 +312,7 @@ struct AngularCorrelationsInJets { if (track.sign() < 0) return isAntiproton; + //TPC if (track.pt() < fAntiprotonTPCTOFpT && TMath::Abs(track.tpcNSigmaPr()) < fAntiprotonTPCnsigLow) isAntiproton = true; if (track.pt() > fAntiprotonTPCTOFpT && TMath::Abs(track.tpcNSigmaPr()) < fAntiprotonTPCnsigHigh) @@ -314,6 +325,7 @@ struct AngularCorrelationsInJets { if (TMath::Abs(track.dcaZ()) > fAntiprotonDCAz) return false; + //TOF if (track.pt() < fAntiprotonTPCTOFpT && TMath::Abs(track.tofNSigmaPr()) < fAntiprotonTOFnsigLow) isAntiproton = true; if (track.pt() > fAntiprotonTPCTOFpT && TMath::Abs(track.tofNSigmaPr()) < fAntiprotonTOFnsigHigh) @@ -329,6 +341,7 @@ struct AngularCorrelationsInJets { if (track.sign() < 0) return isDeuteron; + //TPC if (track.pt() < fDeuteronTPCTOFpT && TMath::Abs(track.tpcNSigmaDe()) < fDeuteronTPCnsigLow) isDeuteron = true; if (track.pt() > fDeuteronTPCTOFpT && TMath::Abs(track.tpcNSigmaDe()) < fDeuteronTPCnsigHigh) @@ -341,6 +354,7 @@ struct AngularCorrelationsInJets { if (TMath::Abs(track.dcaZ()) > fDeuteronDCAz) return false; + //TOF if (track.pt() < fDeuteronTPCTOFpT && TMath::Abs(track.tofNSigmaDe()) < fDeuteronTOFnsigLow) isDeuteron = true; if (track.pt() > fDeuteronTPCTOFpT && TMath::Abs(track.tofNSigmaDe()) < fDeuteronTOFnsigHigh) @@ -356,6 +370,7 @@ struct AngularCorrelationsInJets { if (track.sign() < 0) return isAntideuteron; + //TPC if (track.pt() < fAntideuteronTPCTOFpT && TMath::Abs(track.tpcNSigmaDe()) < fAntideuteronTPCnsigLow) isAntideuteron = true; if (track.pt() > fAntideuteronTPCTOFpT && TMath::Abs(track.tpcNSigmaDe()) < fAntideuteronTPCnsigHigh) @@ -368,6 +383,7 @@ struct AngularCorrelationsInJets { if (TMath::Abs(track.dcaZ()) > fAntideuteronDCAz) return false; + //TOF if (track.pt() < fAntideuteronTPCTOFpT && TMath::Abs(track.tofNSigmaDe()) < fAntideuteronTOFnsigLow) isAntideuteron = true; if (track.pt() > fAntideuteronTPCTOFpT && TMath::Abs(track.tofNSigmaDe()) < fAntideuteronTOFnsigHigh) @@ -376,134 +392,117 @@ struct AngularCorrelationsInJets { return isAntideuteron; } - template - void setTrackBuffer(const T& track, int particleType) //TODO: buffer as argument, lose switch statement + void setTrackBuffer(const auto& tempBuffer, auto& buffer) { - switch (particleType) { - case 1: - if (static_cast(fTrackBufferProton.size()) == fTrackBufferSize) { - fTrackBufferProton.insert(fTrackBufferProton.begin(), track); - fTrackBufferProton.resize(fTrackBufferSize); - } else if (static_cast(fTrackBufferProton.size()) < fTrackBufferSize) { - fTrackBufferProton.emplace_back(track); - } - break; - case 2: - if (static_cast(fTrackBufferAntiproton.size()) == fTrackBufferSize) { - fTrackBufferAntiproton.insert(fTrackBufferAntiproton.begin(), track); - fTrackBufferAntiproton.resize(fTrackBufferSize); - } else if (static_cast(fTrackBufferAntiproton.size()) < fTrackBufferSize) { - fTrackBufferAntiproton.emplace_back(track); - } - break; - case 3: - if (static_cast(fTrackBufferDeuteron.size()) == fTrackBufferSize) { - fTrackBufferDeuteron.insert(fTrackBufferDeuteron.begin(), track); - fTrackBufferDeuteron.resize(fTrackBufferSize); - } else if (static_cast(fTrackBufferDeuteron.size()) < fTrackBufferSize) { - fTrackBufferDeuteron.emplace_back(track); - } - break; - case 4: - if (static_cast(fTrackBufferAntideuteron.size()) == fTrackBufferSize) { - fTrackBufferAntideuteron.insert(fTrackBufferAntideuteron.begin(), track); - fTrackBufferAntideuteron.resize(fTrackBufferSize); - } else if (static_cast(fTrackBufferAntideuteron.size()) < fTrackBufferSize) { - fTrackBufferAntideuteron.emplace_back(track); - } - break; - default: - LOG(warn) << "SetTrackBuffer: invalid particle ID!"; + for (const auto& track : tempBuffer) { + if (static_cast(buffer.size()) == fTrackBufferSize) { + buffer.insert(buffer.begin(), track); + buffer.resize(fTrackBufferSize); + } else if (static_cast(buffer.size()) < fTrackBufferSize) { + buffer.emplace_back(track); + } } } - template - void fillMixedEventDeltas(const T& track, int particleType) //TODO: buffer as argument, lose switch statement + void fillMixedEventDeltas(const auto& track, const auto& buffer, int particleType) { - switch (particleType) { - case 1: - if (fTrackBufferProton.size() == 0) - return; - for (int i = 0; i < static_cast(fTrackBufferProton.size()); i++) { // even if the track buffer isn't even full yet? - double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - fTrackBufferProton[i].phi()); - double DeltaEta = TMath::Abs(track.eta() - fTrackBufferProton[i].eta()); - registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaME"), particleType, DeltaPhi, DeltaEta); - } - break; - case 2: - if (fTrackBufferAntiproton.size() == 0) - return; - for (int i = 0; i < static_cast(fTrackBufferAntiproton.size()); i++) { - double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - fTrackBufferAntiproton[i].phi()); - double DeltaEta = TMath::Abs(track.eta() - fTrackBufferAntiproton[i].eta()); - registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaME"), particleType, DeltaPhi, DeltaEta); - } - break; - case 3: - if (fTrackBufferDeuteron.size() == 0) - return; - for (int i = 0; i < static_cast(fTrackBufferDeuteron.size()); i++) { - double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - fTrackBufferDeuteron[i].phi()); - double DeltaEta = TMath::Abs(track.eta() - fTrackBufferDeuteron[i].eta()); - registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaME"), particleType, DeltaPhi, DeltaEta); - } - break; - case 4: - if (fTrackBufferAntideuteron.size() == 0) - return; - for (int i = 0; i < static_cast(fTrackBufferAntideuteron.size()); i++) { - double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - fTrackBufferAntideuteron[i].phi()); - double DeltaEta = TMath::Abs(track.eta() - fTrackBufferAntideuteron[i].eta()); - registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaME"), particleType, DeltaPhi, DeltaEta); - } - break; - default: - LOG(warn) << "FillMixedEventDeltas: invalid particle ID!"; + if (buffer.size() == 0) + return; + for (int i = 0; i < static_cast(buffer.size()); i++) { // even if the track buffer isn't even full yet? + if (buffer[i].phi() > 2* M_PI || buffer[i].phi() < -2*M_PI) {// problematic? + registryData.fill(HIST("hTrackProtocol"), 11); + continue; + } + double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - buffer[i].phi()); + // LOG(info) << "Lars' Fantastic Debugging Tool | Buffer Track Size/Phi: " << buffer.size() << " | " << i << " | " << buffer[i].phi(); + double DeltaEta = TMath::Abs(track.eta() - buffer[i].eta()); + registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); + registryData.fill(HIST("hDeltaPhiEtaME"), particleType, DeltaPhi, DeltaEta); } } + void doCorrelations(const auto& particleVector, const auto& buffer, auto& tempBuffer, int particleType) + { + for (int i = 0; i < static_cast(particleVector.size()); i++) { + for (int j = i + 1; j < static_cast(particleVector.size()); j++) { + double DeltaPhi = TVector2::Phi_0_2pi(particleVector[i].phi() - particleVector[j].phi()); + double DeltaEta = TMath::Abs(particleVector[i].eta() - particleVector[j].eta()); + if (DeltaPhi > (1.5 * TMath::Pi())) { + DeltaPhi = DeltaPhi - 2 * TMath::Pi(); + } + registryData.fill(HIST("hDeltaPhiSE"), 1, DeltaPhi); + registryData.fill(HIST("hDeltaPhiEtaSE"), 1, DeltaPhi, DeltaEta); + } + fillMixedEventDeltas(particleVector[i], buffer, particleType); + tempBuffer.emplace_back(particleVector[i]); + } + } + template - void fillHistogramsRun2(T const& collision, U const& allTracks) { + void fillHistogramsRun2(T const& collision, U const& allTracks) + { + std::vector fTempBufferProton; + std::vector fTempBufferAntiproton; + std::vector fTempBufferDeuteron; + std::vector fTempBufferAntideuteron; + fTempBufferProton.clear(); + fTempBufferAntiproton.clear(); + fTempBufferDeuteron.clear(); + fTempBufferAntideuteron.clear(); std::vector jetInput; - std::vector constituents; + std::map particles; jetInput.clear(); - constituents.clear(); + particles.clear(); + int index = 0; auto tracks = allTracks.sliceBy(perCollisionFullTracksRun2, collision.globalIndex()); for (const auto& track : tracks) { if (!selectTrack(track)) continue; + + double mass; + if (track.hasTOF()) { + mass = track.mass(); //check reliability, maybe use only pion mass + registryQA.fill(HIST("hTOFmass"), track.pt(), track.mass()); + registryData.fill(HIST("hTrackProtocol"), 4); + } else { + mass = 0.139; //pion mass as default, ~80% are pions + registryData.fill(HIST("hTrackProtocol"), 5); + } registryData.fill(HIST("hPtFullEvent"), track.pt()); - fastjet::PseudoJet inputPseudoJet(track.px(), track.py(), track.pz(), track.energy(track.mass())); - inputPseudoJet.set_user_index(track.globalIndex()+1); + fastjet::PseudoJet inputPseudoJet(track.px(), track.py(), track.pz(), track.energy(mass)); + inputPseudoJet.set_user_index(index); + particles[index] = track; jetInput.emplace_back(inputPseudoJet); + + index++; } if (jetInput.size() < 2) return; registryData.fill(HIST("hEventProtocol"), 2); - fillConstituents(jetInput, constituents); + std::vector constituents = findJets(jetInput); + if (constituents.empty()) + return; - std::vector jetProtons; + std::vector jetProtons; //replace with IDs? std::vector jetAntiprotons; std::vector jetDeuterons; std::vector jetAntideuterons; for (int i = 0; i < static_cast(constituents.size()); i++) { + registryData.fill(HIST("hTrackProtocol"), 10); fastjet::PseudoJet pseudoParticle = constituents[i]; - uint64_t id = pseudoParticle.user_index(); - auto jetParticle = tracks.iteratorAt(id); + int id = pseudoParticle.user_index(); + const auto& jetParticle = particles[id]; + registryData.fill(HIST("hDCAxyFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaXY()); registryData.fill(HIST("hDCAzFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaZ()); registryData.fill(HIST("hTPCsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcSignal()); - registryData.fill(HIST("hTOFsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofSignal()); + registryData.fill(HIST("hTOFsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.beta()); double ptDiff = pseudoParticle.pt() - jetParticle.pt(); registryData.fill(HIST("hPtDiff"), ptDiff); @@ -541,72 +540,29 @@ struct AngularCorrelationsInJets { if ((jetProtons.size() < 2) && (jetAntiprotons.size() < 2) && (jetDeuterons.size() < 2) && (jetAntideuterons.size() < 2)) return; - registryData.fill(HIST("hEventProtocol"), 5); + registryData.fill(HIST("hEventProtocol"), 6); - if (jetProtons.size() > 1) { //TODO: write function that takes jetParticle vector, lose 3/4 if statements - for (int i = 0; i < static_cast(jetProtons.size()); i++) { - for (int j = i + 1; j < static_cast(jetProtons.size()); j++) { - double DeltaPhi = TVector2::Phi_0_2pi(jetProtons[i].phi() - jetProtons[j].phi()); - double DeltaEta = TMath::Abs(jetProtons[i].eta() - jetProtons[j].eta()); - if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); - } - registryData.fill(HIST("hDeltaPhiSE"), 1, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaSE"), 1, DeltaPhi, DeltaEta); - } - fillMixedEventDeltas(jetProtons[i], 1); - setTrackBuffer(jetProtons[i], 1); - } + if (jetProtons.size() > 1) {// unravel functions if code doesn't work? + doCorrelations(jetProtons, fTrackBufferProton, fTempBufferProton, 1); + setTrackBuffer(fTempBufferProton, fTrackBufferProton); } if (jetAntiprotons.size() > 1) { - for (int i = 0; i < static_cast(jetAntiprotons.size()); i++) { - for (int j = i + 1; j < static_cast(jetAntiprotons.size()); j++) { - double DeltaPhi = TVector2::Phi_0_2pi(jetAntiprotons[i].phi() - jetAntiprotons[j].phi()); - double DeltaEta = TMath::Abs(jetAntiprotons[i].eta() - jetAntiprotons[j].eta()); - if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); - } - registryData.fill(HIST("hDeltaPhiSE"), 2, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaSE"), 2, DeltaPhi, DeltaEta); - } - fillMixedEventDeltas(jetAntiprotons[i], 2); - setTrackBuffer(jetAntiprotons[i], 2); - } + doCorrelations(jetAntiprotons, fTrackBufferAntiproton, fTempBufferAntiproton, 2); + setTrackBuffer(fTempBufferAntiproton, fTrackBufferAntiproton); } if (jetDeuterons.size() > 1) { - for (int i = 0; i < static_cast(jetDeuterons.size()); i++) { - for (int j = i + 1; j < static_cast(jetDeuterons.size()); j++) { - double DeltaPhi = TVector2::Phi_0_2pi(jetDeuterons[i].phi() - jetDeuterons[j].phi()); - double DeltaEta = TMath::Abs(jetDeuterons[i].eta() - jetDeuterons[j].eta()); - if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); - } - registryData.fill(HIST("hDeltaPhiSE"), 3, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaSE"), 3, DeltaPhi, DeltaEta); - } - fillMixedEventDeltas(jetDeuterons[i], 3); - setTrackBuffer(jetDeuterons[i], 3); - } + doCorrelations(jetDeuterons, fTrackBufferDeuteron, fTempBufferDeuteron, 3); + setTrackBuffer(fTempBufferDeuteron, fTrackBufferDeuteron); } if (jetAntideuterons.size() > 1) { - for (int i = 0; i < static_cast(jetAntideuterons.size()); i++) { - for (int j = i + 1; j < static_cast(jetAntideuterons.size()); j++) { - double DeltaPhi = TVector2::Phi_0_2pi(jetAntideuterons[i].phi() - jetAntideuterons[j].phi()); - double DeltaEta = TMath::Abs(jetAntideuterons[i].eta() - jetAntideuterons[j].eta()); - if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); - } - registryData.fill(HIST("hDeltaPhiSE"), 4, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaSE"), 4, DeltaPhi, DeltaEta); - } - fillMixedEventDeltas(jetAntideuterons[i], 4); - setTrackBuffer(jetAntideuterons[i], 4); - } + doCorrelations(jetAntideuterons, fTrackBufferAntideuteron, fTempBufferAntideuteron, 4); + setTrackBuffer(fTempBufferAntideuteron, fTrackBufferAntideuteron); } } template - void fillHistogramsRun3(T const& collision, U const& tracks) { + void fillHistogramsRun3(T const& collision, U const& tracks) + { std::vector jetInput; std::vector constituents; jetInput.clear(); @@ -622,8 +578,8 @@ struct AngularCorrelationsInJets { jetInput.emplace_back(inputPseudoJet); } - fillConstituents(jetInput, constituents); - //TODO: add rest of the logic when it is are optimised + findJets(jetInput); + //TODO: add rest of the logic when it is optimised } void processRun2(soa::Join const& collisions, From 75c257031e165aa54de6eb423b53222cb3c24155 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lars=20J=C3=B6rgensen?= Date: Mon, 1 Jul 2024 17:11:22 +0200 Subject: [PATCH 3/8] turned repetitive code into functions (commit to be able to rebase) - there is an issue with this commit --- PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx index 2ccf57883eb..7caa8a81ea7 100644 --- a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx +++ b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx @@ -42,8 +42,6 @@ using namespace o2::framework::expressions; struct AngularCorrelationsInJets { // Preliminary Cuts - Configurable eventSelections{"eventSelection", "sel8", "choose event selection"}; - Configurable trackSelections{"trackSelection", "globalTracks", "set track selection"}; Configurable fMinNCrossedRowsTPC{"minNCrossedRowsTPC", 70, "min number of crossed rows TPC"}; Configurable fMinReqClusterITS{"minReqClusterITS", 2, "min number of clusters required in ITS"}; Configurable fMinReqClusterTPC{"minReqClusterTPC", 70, "min number of clusters required in TPC"}; @@ -392,7 +390,7 @@ struct AngularCorrelationsInJets { return isAntideuteron; } - void setTrackBuffer(const auto& tempBuffer, auto& buffer) + void setTrackBuffer(const auto& tempBuffer, auto& buffer)// tempBuffer problematic? { for (const auto& track : tempBuffer) { if (static_cast(buffer.size()) == fTrackBufferSize) { @@ -409,10 +407,10 @@ struct AngularCorrelationsInJets { if (buffer.size() == 0) return; for (int i = 0; i < static_cast(buffer.size()); i++) { // even if the track buffer isn't even full yet? - if (buffer[i].phi() > 2* M_PI || buffer[i].phi() < -2*M_PI) {// problematic? + /* if (buffer[i].phi() > 2* M_PI || buffer[i].phi() < -2*M_PI) {// problematic? registryData.fill(HIST("hTrackProtocol"), 11); continue; - } + } */ double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - buffer[i].phi()); // LOG(info) << "Lars' Fantastic Debugging Tool | Buffer Track Size/Phi: " << buffer.size() << " | " << i << " | " << buffer[i].phi(); double DeltaEta = TMath::Abs(track.eta() - buffer[i].eta()); @@ -460,7 +458,7 @@ struct AngularCorrelationsInJets { for (const auto& track : tracks) { if (!selectTrack(track)) continue; - + double mass; if (track.hasTOF()) { mass = track.mass(); //check reliability, maybe use only pion mass From baa1d6978197a31a0b33e26de82241c802700724 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lars=20J=C3=B6rgensen?= Date: Wed, 3 Jul 2024 13:27:04 +0200 Subject: [PATCH 4/8] AngularCorrelationsInJets.cxx: add QA histograms --- .../Nuspex/AngularCorrelationsInJets.cxx | 67 ++++++++++++------- 1 file changed, 42 insertions(+), 25 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx index 7caa8a81ea7..9ddb646cbbf 100644 --- a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx +++ b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx @@ -9,7 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -/// \author Lars Jörgensen +// author: Lars Jörgensen #include @@ -50,7 +50,7 @@ struct AngularCorrelationsInJets { Configurable fMaxChi2TPC{"maxChi2TPC", 4.0f, "max chi2 per cluster TPC"}; Configurable fMaxDCAxy{"maxDCA_xy", 0.5f, "max DCA to vertex xy"}; Configurable fMaxDCAz{"maxDCA_z", 2.4f, "max DCA to vertex z"}; - Configurable fMaxEta{"maxEta", 0.8, "max pseudorapidity"}; + Configurable fMaxEta{"maxEta", 0.8, "max pseudorapidity"}; //consider jet cone? // Jet Cuts Configurable fJetR{"jetR", 0.4, "jet resolution parameter"}; @@ -97,12 +97,12 @@ struct AngularCorrelationsInJets { Configurable fTrackBufferSize{"trackBufferSize", 2000, "Number of mixed-event tracks being stored"}; Service ccdb; - // JetFinder jetFinder; - // std::vector constituents; int mRunNumber; - using FullTracksRun2 = soa::Join; - using FullTracksRun3 = soa::Join; + using FullTracksRun2 = soa::Join; + using FullTracksRun3 = soa::Join; using BCsWithRun2Info = soa::Join; Filter prelimTrackCuts = (aod::track::itsChi2NCl < fMaxChi2ITS && @@ -141,11 +141,9 @@ struct AngularCorrelationsInJets { registryData.add("hNumPartInJet", "Number of particles in a jet", HistType::kTH1I, {{200, 0, 200}}); // (Pseudo)Rapidity - registryData.add("hEtaFullEvent", "Particle pseudorapidity;#eta", HistType::kTH1F, {{200, -1, 1}}); registryData.add("hJetRapidity", "Jet rapidity;#it{y}", HistType::kTH1F, {{200, -1, 1}}); // pT - registryData.add("hPtFullEvent", "p_{T} after basic cuts", HistType::kTH1F, {ptAxis}); registryData.add("hPtJetParticle", "p_{T} of particles in jets", HistType::kTH1D, {ptAxis}); registryData.add("hPtSubtractedJet", "Subtracted jet p_{T}", HistType::kTH1D, {ptAxis}); registryData.add("hPtJetProtonDeuteron", "p_{T} of (anti)p, (anti)d", HistType::kTH2D, {particleTypeAxis, ptAxis}); @@ -176,9 +174,19 @@ struct AngularCorrelationsInJets { registryData.add("hJetConeRadius", "Jet Radius;#it{R}", HistType::kTH1F, {{100, 0, 1}}); - //TODO: add QA histograms + // QA registryQA.add("hTOFmass", "TOF mass", HistType::kTH2F, {ptAxis, {1000, 0, 5, "#it{m} [GeV/#it{c}^{2}]"}}); + registryQA.add("hPtFullEvent", "p_{T} after basic cuts", HistType::kTH1F, {ptAxis}); + registryQA.add("hEtaFullEvent", "Particle pseudorapidity;#eta", HistType::kTH1F, {{200, -1, 1}}); registryQA.get(HIST("hTOFmass"))->Sumw2(); + registryQA.add("hCrossedRowsTPC", "Crossed rows TPC", HistType::kTH2I, {ptAxis, {135, 65, 200}}); + registryQA.add("hClusterITS", "ITS clusters", HistType::kTH2I, {ptAxis, {10,0,10}}); + registryQA.add("hClusterTPC", "TPC clusters", HistType::kTH2I, {ptAxis, {135, 65, 200}}); + registryQA.add("hRatioCrossedRowsTPC", "Ratio crossed rows/findable TPC", HistType::kTH2F, {ptAxis, {100,0.5,1.5}}); + registryQA.add("hChi2ITS", "ITS #chi^{2}", HistType::kTH2F, {ptAxis, {400,0,40}}); + registryQA.add("hChi2TPC", "TPC #chi^{2}", HistType::kTH2F, {ptAxis, {50,0,5}}); + registryQA.add("hDCAxyFullEvent", "DCA_{xy} of full event", HistType::kTH2F, {ptAxis, dcaxyAxis}); + registryQA.add("hDCAzFullEvent", "DCA_{z} of full event", HistType::kTH2F, {ptAxis, dcazAxis}); } @@ -232,7 +240,7 @@ struct AngularCorrelationsInJets { fastjet::JetDefinition jetDefBkg(fastjet::kt_algorithm, 0.5); fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(ghost_maxrap, ghost_repeat, ghost_area)); fastjet::AreaDefinition areaDefBkg(fastjet::active_area_explicit_ghosts, fastjet::GhostedAreaSpec(ghost_maxrap)); - fastjet::ClusterSequenceArea clusterSeq(jetInput, jetDef, areaDef); // or CSActiveArea? + fastjet::ClusterSequenceArea clusterSeq(jetInput, jetDef, areaDef); jets = sorted_by_pt(clusterSeq.inclusive_jets()); if (jets.size() == 0) return constituents; @@ -250,16 +258,17 @@ struct AngularCorrelationsInJets { registryData.fill(HIST("hNumberOfJets"), 0); registryData.fill(HIST("hPtTotalJet"), hardestJet.pt()); registryData.fill(HIST("hJetRapidity"), hardestJet.rap()); + registryData.fill(HIST("hNumPartInJet"), hardestJet.constituents().size()); for (const auto& constituent : hardestJet.constituents()) { registryData.fill(HIST("hPtJetParticle"), constituent.pt()); double DeltaPhi = TVector2::Phi_0_2pi(constituent.phi() - hardestJet.phi()); double DeltaEta = constituent.eta() - hardestJet.eta(); - double Delta = TMath::Sqrt(DeltaPhi * DeltaPhi + DeltaEta * DeltaEta); // need 1/pT^2? + double Delta = TMath::Sqrt(DeltaPhi * DeltaPhi + DeltaEta * DeltaEta); registryData.fill(HIST("hJetConeRadius"), Delta); } - fastjet::Selector selector = fastjet::SelectorAbsEtaMax(1.0) * (!fastjet::SelectorNHardest(2)); // fix! + fastjet::Selector selector = fastjet::SelectorAbsEtaMax(1.0) * (!fastjet::SelectorNHardest(2)); //TODO: fix subtraction fastjet::JetMedianBackgroundEstimator bkgEst(selector, jetDefBkg, areaDefBkg); fastjet::Subtractor subtractor(&bkgEst); subtractor.set_use_rho_m(true); @@ -390,7 +399,7 @@ struct AngularCorrelationsInJets { return isAntideuteron; } - void setTrackBuffer(const auto& tempBuffer, auto& buffer)// tempBuffer problematic? + void setTrackBuffer(const auto& tempBuffer, auto& buffer) { for (const auto& track : tempBuffer) { if (static_cast(buffer.size()) == fTrackBufferSize) { @@ -407,12 +416,11 @@ struct AngularCorrelationsInJets { if (buffer.size() == 0) return; for (int i = 0; i < static_cast(buffer.size()); i++) { // even if the track buffer isn't even full yet? - /* if (buffer[i].phi() > 2* M_PI || buffer[i].phi() < -2*M_PI) {// problematic? + if (buffer[i].phi() > 2* M_PI || buffer[i].phi() < -2*M_PI) { registryData.fill(HIST("hTrackProtocol"), 11); continue; - } */ + } double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - buffer[i].phi()); - // LOG(info) << "Lars' Fantastic Debugging Tool | Buffer Track Size/Phi: " << buffer.size() << " | " << i << " | " << buffer[i].phi(); double DeltaEta = TMath::Abs(track.eta() - buffer[i].eta()); registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); registryData.fill(HIST("hDeltaPhiEtaME"), particleType, DeltaPhi, DeltaEta); @@ -469,14 +477,23 @@ struct AngularCorrelationsInJets { registryData.fill(HIST("hTrackProtocol"), 5); } - registryData.fill(HIST("hPtFullEvent"), track.pt()); + registryQA.fill(HIST("hPtFullEvent"), track.pt()); + registryQA.fill(HIST("hEtaFullEvent"), track.eta()); + registryQA.fill(HIST("hCrossedRowsTPC"), track.pt(), track.tpcNClsCrossedRows()); + registryQA.fill(HIST("hClusterITS"), track.pt(), track.itsNCls()); + registryQA.fill(HIST("hClusterTPC"), track.pt(), track.tpcNClsFound()); + registryQA.fill(HIST("hRatioCrossedRowsTPC"), track.pt(), (track.tpcNClsCrossedRows()/track.tpcNClsFindable())); + registryQA.fill(HIST("hChi2ITS"), track.pt(), track.itsChi2NCl()); + registryQA.fill(HIST("hChi2TPC"), track.pt(), track.tpcChi2NCl()); + registryQA.fill(HIST("hDCAxyFullEvent"), track.pt(), track.dcaXY()); + registryQA.fill(HIST("hDCAzFullEvent"), track.pt(), track.dcaZ()); fastjet::PseudoJet inputPseudoJet(track.px(), track.py(), track.pz(), track.energy(mass)); inputPseudoJet.set_user_index(index); particles[index] = track; jetInput.emplace_back(inputPseudoJet); index++; - } + } //for (const auto& track : tracks) if (jetInput.size() < 2) return; @@ -500,7 +517,8 @@ struct AngularCorrelationsInJets { registryData.fill(HIST("hDCAxyFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaXY()); registryData.fill(HIST("hDCAzFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaZ()); registryData.fill(HIST("hTPCsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcSignal()); - registryData.fill(HIST("hTOFsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.beta()); + if (jetParticle.hasTOF()) + registryData.fill(HIST("hTOFsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.beta()); double ptDiff = pseudoParticle.pt() - jetParticle.pt(); registryData.fill(HIST("hPtDiff"), ptDiff); @@ -510,7 +528,8 @@ struct AngularCorrelationsInJets { continue; if (isProton(jetParticle) || isAntiproton(jetParticle)) { // collect (anti)protons in jet registryData.fill(HIST("hTPCnsigmaProton"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcNSigmaPr()); - registryData.fill(HIST("hTOFnsigmaProton"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaPr()); + if (jetParticle.hasTOF()) + registryData.fill(HIST("hTOFnsigmaProton"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaPr()); if (isProton(jetParticle)) { particleType = 1; registryData.fill(HIST("hTrackProtocol"), 6); // # protons @@ -522,7 +541,8 @@ struct AngularCorrelationsInJets { } } else if (isDeuteron(jetParticle) || isAntideuteron(jetParticle)) { // collect (anti)deuterons in jet registryData.fill(HIST("hTPCnsigmaDeuteron"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcNSigmaDe()); - registryData.fill(HIST("hTOFnsigmaDeuteron"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaDe()); + if (jetParticle.hasTOF()) + registryData.fill(HIST("hTOFnsigmaDeuteron"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaDe()); if (isDeuteron(jetParticle)) { particleType = 3; registryData.fill(HIST("hTrackProtocol"), 8); // # deuterons @@ -540,7 +560,7 @@ struct AngularCorrelationsInJets { return; registryData.fill(HIST("hEventProtocol"), 6); - if (jetProtons.size() > 1) {// unravel functions if code doesn't work? + if (jetProtons.size() > 1) { doCorrelations(jetProtons, fTrackBufferProton, fTempBufferProton, 1); setTrackBuffer(fTempBufferProton, fTrackBufferProton); } @@ -576,7 +596,6 @@ struct AngularCorrelationsInJets { jetInput.emplace_back(inputPseudoJet); } - findJets(jetInput); //TODO: add rest of the logic when it is optimised } @@ -584,9 +603,7 @@ struct AngularCorrelationsInJets { soa::Filtered const& tracks, BCsWithRun2Info const&) { - int index = 0; for (const auto& collision : collisions) { - index++; auto bc = collision.bc_as(); initCCDB(bc); From a764b49bfc758a7450bb02fc60c6d5fe0e12059f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lars=20J=C3=B6rgensen?= Date: Wed, 10 Jul 2024 13:16:27 +0200 Subject: [PATCH 5/8] AngularCorrelationsInJets.cxx: fixed issues and added Run3 logic --- .../Nuspex/AngularCorrelationsInJets.cxx | 208 +++++++++++++++--- 1 file changed, 173 insertions(+), 35 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx index 9ddb646cbbf..f905ad9fb15 100644 --- a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx +++ b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx @@ -102,7 +102,7 @@ struct AngularCorrelationsInJets { using FullTracksRun2 = soa::Join; using FullTracksRun3 = soa::Join; + aod::TracksDCA, aod::pidTPCFullPr, aod::pidTPCFullDe, aod::pidTOFFullPr, aod::pidTOFFullDe, aod::pidTOFmass, aod::pidTOFbeta>; using BCsWithRun2Info = soa::Join; Filter prelimTrackCuts = (aod::track::itsChi2NCl < fMaxChi2ITS && @@ -112,6 +112,7 @@ struct AngularCorrelationsInJets { nabs(aod::track::eta) < fMaxEta); Preslice perCollisionFullTracksRun2 = o2::aod::track::collisionId; + Preslice perCollisionFullTracksRun3 = o2::aod::track::collisionId; AxisSpec ptAxis = {1000, 0, 100, "#it{p}_{T} [GeV/#it{c}]"}; AxisSpec particleTypeAxis = {4, 1, 5, "[p, ap, d, ad]"}; @@ -190,10 +191,15 @@ struct AngularCorrelationsInJets { } - std::vector fTrackBufferProton; - std::vector fTrackBufferAntiproton; - std::vector fTrackBufferDeuteron; - std::vector fTrackBufferAntideuteron; + + std::vector fTrackBufferProtonRun2; + std::vector fTrackBufferAntiprotonRun2; + std::vector fTrackBufferDeuteronRun2; + std::vector fTrackBufferAntideuteronRun2; + std::vector fTrackBufferProtonRun3; + std::vector fTrackBufferAntiprotonRun3; + std::vector fTrackBufferDeuteronRun3; + std::vector fTrackBufferAntideuteronRun3; //TODO: check if FullTracksRun2 works for Run3 too or add Run3 track buffers template @@ -415,13 +421,19 @@ struct AngularCorrelationsInJets { { if (buffer.size() == 0) return; - for (int i = 0; i < static_cast(buffer.size()); i++) { // even if the track buffer isn't even full yet? - if (buffer[i].phi() > 2* M_PI || buffer[i].phi() < -2*M_PI) { - registryData.fill(HIST("hTrackProtocol"), 11); + for (int i = 0; i < static_cast(buffer.size()); i++) { + if (std::isnan(buffer.at(i).phi())) continue; + if (buffer.at(i).phi() > 2*TMath::Pi() || buffer.at(i).phi() < -2*TMath::Pi()) { + registryData.fill(HIST("hTrackProtocol"), 13); + continue; + } + + double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - buffer.at(i).phi()); + if (DeltaPhi > (1.5 * TMath::Pi())) { + DeltaPhi = DeltaPhi - 2 * TMath::Pi(); } - double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - buffer[i].phi()); - double DeltaEta = TMath::Abs(track.eta() - buffer[i].eta()); + double DeltaEta = TMath::Abs(track.eta() - buffer.at(i).eta()); registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); registryData.fill(HIST("hDeltaPhiEtaME"), particleType, DeltaPhi, DeltaEta); } @@ -429,18 +441,29 @@ struct AngularCorrelationsInJets { void doCorrelations(const auto& particleVector, const auto& buffer, auto& tempBuffer, int particleType) { - for (int i = 0; i < static_cast(particleVector.size()); i++) { + for (int i = 0; i < static_cast(particleVector.size()); i++) {// maybe simply introduce phi cut? + if (particleVector.at(i).phi() > 2*TMath::Pi() || particleVector.at(i).phi() < -2*TMath::Pi()) { + registryData.fill(HIST("hTrackProtocol"), 11); + continue; + } for (int j = i + 1; j < static_cast(particleVector.size()); j++) { - double DeltaPhi = TVector2::Phi_0_2pi(particleVector[i].phi() - particleVector[j].phi()); - double DeltaEta = TMath::Abs(particleVector[i].eta() - particleVector[j].eta()); + if ((j == static_cast(particleVector.size())) || std::isnan(particleVector.at(j).phi())) + continue; + if (particleVector.at(j).phi() > 2*TMath::Pi() || particleVector.at(j).phi() < -2*TMath::Pi()) { + registryData.fill(HIST("hTrackProtocol"), 12); + continue; + } + + double DeltaPhi = TVector2::Phi_0_2pi(particleVector.at(i).phi() - particleVector.at(j).phi()); + double DeltaEta = TMath::Abs(particleVector.at(i).eta() - particleVector[j].eta()); if (DeltaPhi > (1.5 * TMath::Pi())) { DeltaPhi = DeltaPhi - 2 * TMath::Pi(); } registryData.fill(HIST("hDeltaPhiSE"), 1, DeltaPhi); registryData.fill(HIST("hDeltaPhiEtaSE"), 1, DeltaPhi, DeltaEta); } - fillMixedEventDeltas(particleVector[i], buffer, particleType); - tempBuffer.emplace_back(particleVector[i]); + fillMixedEventDeltas(particleVector.at(i), buffer, particleType); + tempBuffer.emplace_back(particleVector.at(i)); } } @@ -477,12 +500,13 @@ struct AngularCorrelationsInJets { registryData.fill(HIST("hTrackProtocol"), 5); } + // double ratioCrossedRowsTPC = track.tpcNClsCrossedRows()/track.tpcNClsFindable(); registryQA.fill(HIST("hPtFullEvent"), track.pt()); registryQA.fill(HIST("hEtaFullEvent"), track.eta()); registryQA.fill(HIST("hCrossedRowsTPC"), track.pt(), track.tpcNClsCrossedRows()); registryQA.fill(HIST("hClusterITS"), track.pt(), track.itsNCls()); registryQA.fill(HIST("hClusterTPC"), track.pt(), track.tpcNClsFound()); - registryQA.fill(HIST("hRatioCrossedRowsTPC"), track.pt(), (track.tpcNClsCrossedRows()/track.tpcNClsFindable())); + // registryQA.fill(HIST("hRatioCrossedRowsTPC"), track.pt(), ratioCrossedRowsTPC); registryQA.fill(HIST("hChi2ITS"), track.pt(), track.itsChi2NCl()); registryQA.fill(HIST("hChi2TPC"), track.pt(), track.tpcChi2NCl()); registryQA.fill(HIST("hDCAxyFullEvent"), track.pt(), track.dcaXY()); @@ -510,7 +534,7 @@ struct AngularCorrelationsInJets { for (int i = 0; i < static_cast(constituents.size()); i++) { registryData.fill(HIST("hTrackProtocol"), 10); - fastjet::PseudoJet pseudoParticle = constituents[i]; + fastjet::PseudoJet pseudoParticle = constituents.at(i); int id = pseudoParticle.user_index(); const auto& jetParticle = particles[id]; @@ -561,42 +585,156 @@ struct AngularCorrelationsInJets { registryData.fill(HIST("hEventProtocol"), 6); if (jetProtons.size() > 1) { - doCorrelations(jetProtons, fTrackBufferProton, fTempBufferProton, 1); - setTrackBuffer(fTempBufferProton, fTrackBufferProton); + doCorrelations(jetProtons, fTrackBufferProtonRun2, fTempBufferProton, 1); + setTrackBuffer(fTempBufferProton, fTrackBufferProtonRun2); } if (jetAntiprotons.size() > 1) { - doCorrelations(jetAntiprotons, fTrackBufferAntiproton, fTempBufferAntiproton, 2); - setTrackBuffer(fTempBufferAntiproton, fTrackBufferAntiproton); + doCorrelations(jetAntiprotons, fTrackBufferAntiprotonRun2, fTempBufferAntiproton, 2); + setTrackBuffer(fTempBufferAntiproton, fTrackBufferAntiprotonRun2); } if (jetDeuterons.size() > 1) { - doCorrelations(jetDeuterons, fTrackBufferDeuteron, fTempBufferDeuteron, 3); - setTrackBuffer(fTempBufferDeuteron, fTrackBufferDeuteron); + doCorrelations(jetDeuterons, fTrackBufferDeuteronRun2, fTempBufferDeuteron, 3); + setTrackBuffer(fTempBufferDeuteron, fTrackBufferDeuteronRun2); } if (jetAntideuterons.size() > 1) { - doCorrelations(jetAntideuterons, fTrackBufferAntideuteron, fTempBufferAntideuteron, 4); - setTrackBuffer(fTempBufferAntideuteron, fTrackBufferAntideuteron); + doCorrelations(jetAntideuterons, fTrackBufferAntideuteronRun2, fTempBufferAntideuteron, 4); + setTrackBuffer(fTempBufferAntideuteron, fTrackBufferAntideuteronRun2); } } template - void fillHistogramsRun3(T const& collision, U const& tracks) + void fillHistogramsRun3(T const& collision, U const& allTracks) { + std::vector fTempBufferProton; + std::vector fTempBufferAntiproton; + std::vector fTempBufferDeuteron; + std::vector fTempBufferAntideuteron; + fTempBufferProton.clear(); + fTempBufferAntiproton.clear(); + fTempBufferDeuteron.clear(); + fTempBufferAntideuteron.clear(); std::vector jetInput; - std::vector constituents; + std::map particles; jetInput.clear(); - constituents.clear(); + particles.clear(); + int index = 0; + + auto tracks = allTracks.sliceBy(perCollisionFullTracksRun2, collision.globalIndex()); for (const auto& track : tracks) { - registryData.fill(HIST("hTrackProtocol"), 0); - registryData.fill(HIST("hEtaFullEvent"), track.eta()); - registryData.fill(HIST("hPtFullEvent"), track.pt()); + if (!selectTrack(track)) + continue; - fastjet::PseudoJet inputPseudoJet(track.px(), track.py(), track.pz(), track.energy(track.mass())); - inputPseudoJet.set_user_index(track.globalIndex()); + double mass; + if (track.hasTOF()) { + mass = track.mass(); //check reliability, maybe use only pion mass + registryQA.fill(HIST("hTOFmass"), track.pt(), track.mass()); + registryData.fill(HIST("hTrackProtocol"), 4); + } else { + mass = 0.139; //pion mass as default, ~80% are pions + registryData.fill(HIST("hTrackProtocol"), 5); + } + + // double ratioCrossedRowsTPC = track.tpcNClsCrossedRows()/track.tpcNClsFindable(); + registryQA.fill(HIST("hPtFullEvent"), track.pt()); + registryQA.fill(HIST("hEtaFullEvent"), track.eta()); + registryQA.fill(HIST("hCrossedRowsTPC"), track.pt(), track.tpcNClsCrossedRows()); + registryQA.fill(HIST("hClusterITS"), track.pt(), track.itsNCls()); + registryQA.fill(HIST("hClusterTPC"), track.pt(), track.tpcNClsFound()); + // registryQA.fill(HIST("hRatioCrossedRowsTPC"), track.pt(), ratioCrossedRowsTPC); + registryQA.fill(HIST("hChi2ITS"), track.pt(), track.itsChi2NCl()); + registryQA.fill(HIST("hChi2TPC"), track.pt(), track.tpcChi2NCl()); + registryQA.fill(HIST("hDCAxyFullEvent"), track.pt(), track.dcaXY()); + registryQA.fill(HIST("hDCAzFullEvent"), track.pt(), track.dcaZ()); + fastjet::PseudoJet inputPseudoJet(track.px(), track.py(), track.pz(), track.energy(mass)); + inputPseudoJet.set_user_index(index); + particles[index] = track; jetInput.emplace_back(inputPseudoJet); - } - //TODO: add rest of the logic when it is optimised + index++; + } //for (const auto& track : tracks) + + if (jetInput.size() < 2) + return; + registryData.fill(HIST("hEventProtocol"), 2); + + std::vector constituents = findJets(jetInput); + if (constituents.empty()) + return; + + std::vector jetProtons; //replace with IDs? + std::vector jetAntiprotons; + std::vector jetDeuterons; + std::vector jetAntideuterons; + + for (int i = 0; i < static_cast(constituents.size()); i++) { + registryData.fill(HIST("hTrackProtocol"), 10); + fastjet::PseudoJet pseudoParticle = constituents.at(i); + int id = pseudoParticle.user_index(); + const auto& jetParticle = particles[id]; + + registryData.fill(HIST("hDCAxyFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaXY()); + registryData.fill(HIST("hDCAzFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaZ()); + registryData.fill(HIST("hTPCsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcSignal()); + if (jetParticle.hasTOF()) + registryData.fill(HIST("hTOFsignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.beta()); + + double ptDiff = pseudoParticle.pt() - jetParticle.pt(); + registryData.fill(HIST("hPtDiff"), ptDiff); + int particleType = 0; + + if (jetParticle.pt() < fMinJetParticlePt) + continue; + if (isProton(jetParticle) || isAntiproton(jetParticle)) { // collect (anti)protons in jet + registryData.fill(HIST("hTPCnsigmaProton"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcNSigmaPr()); + if (jetParticle.hasTOF()) + registryData.fill(HIST("hTOFnsigmaProton"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaPr()); + if (isProton(jetParticle)) { + particleType = 1; + registryData.fill(HIST("hTrackProtocol"), 6); // # protons + jetProtons.emplace_back(jetParticle); + } else { + particleType = 2; + registryData.fill(HIST("hTrackProtocol"), 7); // # antiprotons + jetAntiprotons.emplace_back(jetParticle); + } + } else if (isDeuteron(jetParticle) || isAntideuteron(jetParticle)) { // collect (anti)deuterons in jet + registryData.fill(HIST("hTPCnsigmaDeuteron"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcNSigmaDe()); + if (jetParticle.hasTOF()) + registryData.fill(HIST("hTOFnsigmaDeuteron"), jetParticle.pt() * jetParticle.sign(), jetParticle.tofNSigmaDe()); + if (isDeuteron(jetParticle)) { + particleType = 3; + registryData.fill(HIST("hTrackProtocol"), 8); // # deuterons + jetDeuterons.emplace_back(jetParticle); + } else { + particleType = 4; + registryData.fill(HIST("hTrackProtocol"), 9); // # antideuterons + jetAntideuterons.emplace_back(jetParticle); + } + } + registryData.fill(HIST("hPtJetProtonDeuteron"), particleType, jetParticle.pt()); + } // for (int i=0; i(constituents.size()); i++) + + if ((jetProtons.size() < 2) && (jetAntiprotons.size() < 2) && (jetDeuterons.size() < 2) && (jetAntideuterons.size() < 2)) + return; + registryData.fill(HIST("hEventProtocol"), 6); + + if (jetProtons.size() > 1) { + doCorrelations(jetProtons, fTrackBufferProtonRun3, fTempBufferProton, 1); + setTrackBuffer(fTempBufferProton, fTrackBufferProtonRun3); + } + if (jetAntiprotons.size() > 1) { + doCorrelations(jetAntiprotons, fTrackBufferAntiprotonRun3, fTempBufferAntiproton, 2); + setTrackBuffer(fTempBufferAntiproton, fTrackBufferAntiprotonRun3); + } + if (jetDeuterons.size() > 1) { + doCorrelations(jetDeuterons, fTrackBufferDeuteronRun3, fTempBufferDeuteron, 3); + setTrackBuffer(fTempBufferDeuteron, fTrackBufferDeuteronRun3); + } + if (jetAntideuterons.size() > 1) { + doCorrelations(jetAntideuterons, fTrackBufferAntideuteronRun3, fTempBufferAntideuteron, 4); + setTrackBuffer(fTempBufferAntideuteron, fTrackBufferAntideuteronRun3); + } } void processRun2(soa::Join const& collisions, From 5ba1d8f2e739fdd1eed91ed5750e3b20849fd836 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 10 Jul 2024 11:19:19 +0000 Subject: [PATCH 6/8] Please consider the following formatting changes --- .../Nuspex/AngularCorrelationsInJets.cxx | 103 +++++++++--------- 1 file changed, 50 insertions(+), 53 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx index f905ad9fb15..a74564dde81 100644 --- a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx +++ b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx @@ -50,7 +50,7 @@ struct AngularCorrelationsInJets { Configurable fMaxChi2TPC{"maxChi2TPC", 4.0f, "max chi2 per cluster TPC"}; Configurable fMaxDCAxy{"maxDCA_xy", 0.5f, "max DCA to vertex xy"}; Configurable fMaxDCAz{"maxDCA_z", 2.4f, "max DCA to vertex z"}; - Configurable fMaxEta{"maxEta", 0.8, "max pseudorapidity"}; //consider jet cone? + Configurable fMaxEta{"maxEta", 0.8, "max pseudorapidity"}; // consider jet cone? // Jet Cuts Configurable fJetR{"jetR", 0.4, "jet resolution parameter"}; @@ -100,7 +100,7 @@ struct AngularCorrelationsInJets { int mRunNumber; using FullTracksRun2 = soa::Join; + aod::TrackSelectionExtension, aod::TracksDCA, aod::pidTPCFullPr, aod::pidTPCFullDe, aod::pidTOFFullPr, aod::pidTOFFullDe, aod::pidTOFmass, aod::pidTOFbeta>; using FullTracksRun3 = soa::Join; using BCsWithRun2Info = soa::Join; @@ -128,7 +128,7 @@ struct AngularCorrelationsInJets { void init(o2::framework::InitContext&) { - mRunNumber = 0; + mRunNumber = 0; ccdb->setURL("http://alice-ccdb.cern.ch"); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -181,17 +181,15 @@ struct AngularCorrelationsInJets { registryQA.add("hEtaFullEvent", "Particle pseudorapidity;#eta", HistType::kTH1F, {{200, -1, 1}}); registryQA.get(HIST("hTOFmass"))->Sumw2(); registryQA.add("hCrossedRowsTPC", "Crossed rows TPC", HistType::kTH2I, {ptAxis, {135, 65, 200}}); - registryQA.add("hClusterITS", "ITS clusters", HistType::kTH2I, {ptAxis, {10,0,10}}); + registryQA.add("hClusterITS", "ITS clusters", HistType::kTH2I, {ptAxis, {10, 0, 10}}); registryQA.add("hClusterTPC", "TPC clusters", HistType::kTH2I, {ptAxis, {135, 65, 200}}); - registryQA.add("hRatioCrossedRowsTPC", "Ratio crossed rows/findable TPC", HistType::kTH2F, {ptAxis, {100,0.5,1.5}}); - registryQA.add("hChi2ITS", "ITS #chi^{2}", HistType::kTH2F, {ptAxis, {400,0,40}}); - registryQA.add("hChi2TPC", "TPC #chi^{2}", HistType::kTH2F, {ptAxis, {50,0,5}}); + registryQA.add("hRatioCrossedRowsTPC", "Ratio crossed rows/findable TPC", HistType::kTH2F, {ptAxis, {100, 0.5, 1.5}}); + registryQA.add("hChi2ITS", "ITS #chi^{2}", HistType::kTH2F, {ptAxis, {400, 0, 40}}); + registryQA.add("hChi2TPC", "TPC #chi^{2}", HistType::kTH2F, {ptAxis, {50, 0, 5}}); registryQA.add("hDCAxyFullEvent", "DCA_{xy} of full event", HistType::kTH2F, {ptAxis, dcaxyAxis}); registryQA.add("hDCAzFullEvent", "DCA_{z} of full event", HistType::kTH2F, {ptAxis, dcazAxis}); - } - std::vector fTrackBufferProtonRun2; std::vector fTrackBufferAntiprotonRun2; std::vector fTrackBufferDeuteronRun2; @@ -200,7 +198,7 @@ struct AngularCorrelationsInJets { std::vector fTrackBufferAntiprotonRun3; std::vector fTrackBufferDeuteronRun3; std::vector fTrackBufferAntideuteronRun3; - //TODO: check if FullTracksRun2 works for Run3 too or add Run3 track buffers + // TODO: check if FullTracksRun2 works for Run3 too or add Run3 track buffers template void initCCDB(Bc const& bc) @@ -221,9 +219,10 @@ struct AngularCorrelationsInJets { return false; } if (doprocessRun2) { - if (!(track.trackType() & o2::aod::track::Run2Track)/* || - !(track.flags() & o2::aod::track::TPCrefit) || - !(track.flags() & o2::aod::track::ITSrefit) */) { + if (!(track.trackType() & o2::aod::track::Run2Track) /* || + !(track.flags() & o2::aod::track::TPCrefit) || + !(track.flags() & o2::aod::track::ITSrefit) */ + ) { return false; } } @@ -274,7 +273,7 @@ struct AngularCorrelationsInJets { registryData.fill(HIST("hJetConeRadius"), Delta); } - fastjet::Selector selector = fastjet::SelectorAbsEtaMax(1.0) * (!fastjet::SelectorNHardest(2)); //TODO: fix subtraction + fastjet::Selector selector = fastjet::SelectorAbsEtaMax(1.0) * (!fastjet::SelectorNHardest(2)); // TODO: fix subtraction fastjet::JetMedianBackgroundEstimator bkgEst(selector, jetDefBkg, areaDefBkg); fastjet::Subtractor subtractor(&bkgEst); subtractor.set_use_rho_m(true); @@ -296,7 +295,7 @@ struct AngularCorrelationsInJets { if (track.sign() < 0) return isProton; - //TPC + // TPC if (track.pt() < fProtonTPCTOFpT && TMath::Abs(track.tpcNSigmaPr()) < fProtonTPCnsigLow) isProton = true; if (track.pt() > fProtonTPCTOFpT && TMath::Abs(track.tpcNSigmaPr()) < fProtonTPCnsigHigh) @@ -309,7 +308,7 @@ struct AngularCorrelationsInJets { if (TMath::Abs(track.dcaZ()) > fProtonDCAz) return false; - //TOF + // TOF if (track.pt() < fProtonTPCTOFpT && TMath::Abs(track.tofNSigmaPr()) < fProtonTOFnsigLow) isProton = true; if (track.pt() > fProtonTPCTOFpT && TMath::Abs(track.tofNSigmaPr()) < fProtonTOFnsigHigh) @@ -325,7 +324,7 @@ struct AngularCorrelationsInJets { if (track.sign() < 0) return isAntiproton; - //TPC + // TPC if (track.pt() < fAntiprotonTPCTOFpT && TMath::Abs(track.tpcNSigmaPr()) < fAntiprotonTPCnsigLow) isAntiproton = true; if (track.pt() > fAntiprotonTPCTOFpT && TMath::Abs(track.tpcNSigmaPr()) < fAntiprotonTPCnsigHigh) @@ -338,7 +337,7 @@ struct AngularCorrelationsInJets { if (TMath::Abs(track.dcaZ()) > fAntiprotonDCAz) return false; - //TOF + // TOF if (track.pt() < fAntiprotonTPCTOFpT && TMath::Abs(track.tofNSigmaPr()) < fAntiprotonTOFnsigLow) isAntiproton = true; if (track.pt() > fAntiprotonTPCTOFpT && TMath::Abs(track.tofNSigmaPr()) < fAntiprotonTOFnsigHigh) @@ -354,7 +353,7 @@ struct AngularCorrelationsInJets { if (track.sign() < 0) return isDeuteron; - //TPC + // TPC if (track.pt() < fDeuteronTPCTOFpT && TMath::Abs(track.tpcNSigmaDe()) < fDeuteronTPCnsigLow) isDeuteron = true; if (track.pt() > fDeuteronTPCTOFpT && TMath::Abs(track.tpcNSigmaDe()) < fDeuteronTPCnsigHigh) @@ -367,7 +366,7 @@ struct AngularCorrelationsInJets { if (TMath::Abs(track.dcaZ()) > fDeuteronDCAz) return false; - //TOF + // TOF if (track.pt() < fDeuteronTPCTOFpT && TMath::Abs(track.tofNSigmaDe()) < fDeuteronTOFnsigLow) isDeuteron = true; if (track.pt() > fDeuteronTPCTOFpT && TMath::Abs(track.tofNSigmaDe()) < fDeuteronTOFnsigHigh) @@ -383,7 +382,7 @@ struct AngularCorrelationsInJets { if (track.sign() < 0) return isAntideuteron; - //TPC + // TPC if (track.pt() < fAntideuteronTPCTOFpT && TMath::Abs(track.tpcNSigmaDe()) < fAntideuteronTPCnsigLow) isAntideuteron = true; if (track.pt() > fAntideuteronTPCTOFpT && TMath::Abs(track.tpcNSigmaDe()) < fAntideuteronTPCnsigHigh) @@ -396,7 +395,7 @@ struct AngularCorrelationsInJets { if (TMath::Abs(track.dcaZ()) > fAntideuteronDCAz) return false; - //TOF + // TOF if (track.pt() < fAntideuteronTPCTOFpT && TMath::Abs(track.tofNSigmaDe()) < fAntideuteronTOFnsigLow) isAntideuteron = true; if (track.pt() > fAntideuteronTPCTOFpT && TMath::Abs(track.tofNSigmaDe()) < fAntideuteronTOFnsigHigh) @@ -424,14 +423,14 @@ struct AngularCorrelationsInJets { for (int i = 0; i < static_cast(buffer.size()); i++) { if (std::isnan(buffer.at(i).phi())) continue; - if (buffer.at(i).phi() > 2*TMath::Pi() || buffer.at(i).phi() < -2*TMath::Pi()) { + if (buffer.at(i).phi() > 2 * TMath::Pi() || buffer.at(i).phi() < -2 * TMath::Pi()) { registryData.fill(HIST("hTrackProtocol"), 13); continue; } double DeltaPhi = TVector2::Phi_0_2pi(track.phi() - buffer.at(i).phi()); if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); + DeltaPhi = DeltaPhi - 2 * TMath::Pi(); } double DeltaEta = TMath::Abs(track.eta() - buffer.at(i).eta()); registryData.fill(HIST("hDeltaPhiME"), particleType, DeltaPhi); @@ -441,30 +440,30 @@ struct AngularCorrelationsInJets { void doCorrelations(const auto& particleVector, const auto& buffer, auto& tempBuffer, int particleType) { - for (int i = 0; i < static_cast(particleVector.size()); i++) {// maybe simply introduce phi cut? - if (particleVector.at(i).phi() > 2*TMath::Pi() || particleVector.at(i).phi() < -2*TMath::Pi()) { + for (int i = 0; i < static_cast(particleVector.size()); i++) { // maybe simply introduce phi cut? + if (particleVector.at(i).phi() > 2 * TMath::Pi() || particleVector.at(i).phi() < -2 * TMath::Pi()) { registryData.fill(HIST("hTrackProtocol"), 11); continue; } - for (int j = i + 1; j < static_cast(particleVector.size()); j++) { - if ((j == static_cast(particleVector.size())) || std::isnan(particleVector.at(j).phi())) - continue; - if (particleVector.at(j).phi() > 2*TMath::Pi() || particleVector.at(j).phi() < -2*TMath::Pi()) { - registryData.fill(HIST("hTrackProtocol"), 12); - continue; - } - - double DeltaPhi = TVector2::Phi_0_2pi(particleVector.at(i).phi() - particleVector.at(j).phi()); - double DeltaEta = TMath::Abs(particleVector.at(i).eta() - particleVector[j].eta()); - if (DeltaPhi > (1.5 * TMath::Pi())) { - DeltaPhi = DeltaPhi - 2 * TMath::Pi(); - } - registryData.fill(HIST("hDeltaPhiSE"), 1, DeltaPhi); - registryData.fill(HIST("hDeltaPhiEtaSE"), 1, DeltaPhi, DeltaEta); + for (int j = i + 1; j < static_cast(particleVector.size()); j++) { + if ((j == static_cast(particleVector.size())) || std::isnan(particleVector.at(j).phi())) + continue; + if (particleVector.at(j).phi() > 2 * TMath::Pi() || particleVector.at(j).phi() < -2 * TMath::Pi()) { + registryData.fill(HIST("hTrackProtocol"), 12); + continue; + } + + double DeltaPhi = TVector2::Phi_0_2pi(particleVector.at(i).phi() - particleVector.at(j).phi()); + double DeltaEta = TMath::Abs(particleVector.at(i).eta() - particleVector[j].eta()); + if (DeltaPhi > (1.5 * TMath::Pi())) { + DeltaPhi = DeltaPhi - 2 * TMath::Pi(); } - fillMixedEventDeltas(particleVector.at(i), buffer, particleType); - tempBuffer.emplace_back(particleVector.at(i)); + registryData.fill(HIST("hDeltaPhiSE"), 1, DeltaPhi); + registryData.fill(HIST("hDeltaPhiEtaSE"), 1, DeltaPhi, DeltaEta); } + fillMixedEventDeltas(particleVector.at(i), buffer, particleType); + tempBuffer.emplace_back(particleVector.at(i)); + } } template @@ -492,11 +491,11 @@ struct AngularCorrelationsInJets { double mass; if (track.hasTOF()) { - mass = track.mass(); //check reliability, maybe use only pion mass + mass = track.mass(); // check reliability, maybe use only pion mass registryQA.fill(HIST("hTOFmass"), track.pt(), track.mass()); registryData.fill(HIST("hTrackProtocol"), 4); } else { - mass = 0.139; //pion mass as default, ~80% are pions + mass = 0.139; // pion mass as default, ~80% are pions registryData.fill(HIST("hTrackProtocol"), 5); } @@ -517,7 +516,7 @@ struct AngularCorrelationsInJets { jetInput.emplace_back(inputPseudoJet); index++; - } //for (const auto& track : tracks) + } // for (const auto& track : tracks) if (jetInput.size() < 2) return; @@ -527,7 +526,7 @@ struct AngularCorrelationsInJets { if (constituents.empty()) return; - std::vector jetProtons; //replace with IDs? + std::vector jetProtons; // replace with IDs? std::vector jetAntiprotons; std::vector jetDeuterons; std::vector jetAntideuterons; @@ -627,11 +626,11 @@ struct AngularCorrelationsInJets { double mass; if (track.hasTOF()) { - mass = track.mass(); //check reliability, maybe use only pion mass + mass = track.mass(); // check reliability, maybe use only pion mass registryQA.fill(HIST("hTOFmass"), track.pt(), track.mass()); registryData.fill(HIST("hTrackProtocol"), 4); } else { - mass = 0.139; //pion mass as default, ~80% are pions + mass = 0.139; // pion mass as default, ~80% are pions registryData.fill(HIST("hTrackProtocol"), 5); } @@ -652,7 +651,7 @@ struct AngularCorrelationsInJets { jetInput.emplace_back(inputPseudoJet); index++; - } //for (const auto& track : tracks) + } // for (const auto& track : tracks) if (jetInput.size() < 2) return; @@ -662,7 +661,7 @@ struct AngularCorrelationsInJets { if (constituents.empty()) return; - std::vector jetProtons; //replace with IDs? + std::vector jetProtons; // replace with IDs? std::vector jetAntiprotons; std::vector jetDeuterons; std::vector jetAntideuterons; @@ -763,8 +762,6 @@ struct AngularCorrelationsInJets { PROCESS_SWITCH(AngularCorrelationsInJets, processRun3, "process Run 3 data", false); }; - - WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ From 7db978addc7e03122d2bfaee09425dabba05503d Mon Sep 17 00:00:00 2001 From: Lars <146946151+ljoergen@users.noreply.github.com> Date: Wed, 10 Jul 2024 13:31:06 +0200 Subject: [PATCH 7/8] AngularCorrelationsInJets.cxx: changed comment from /**/ to // to stop megalinter from complaining --- PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx index a74564dde81..7a45fdfa525 100644 --- a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx +++ b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx @@ -219,9 +219,9 @@ struct AngularCorrelationsInJets { return false; } if (doprocessRun2) { - if (!(track.trackType() & o2::aod::track::Run2Track) /* || - !(track.flags() & o2::aod::track::TPCrefit) || - !(track.flags() & o2::aod::track::ITSrefit) */ + if (!(track.trackType() & o2::aod::track::Run2Track) //|| + //!(track.flags() & o2::aod::track::TPCrefit) || + //!(track.flags() & o2::aod::track::ITSrefit) ) { return false; } From e277eca838aabf88744f2d7dfcb3796102d6f823 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 10 Jul 2024 11:32:11 +0000 Subject: [PATCH 8/8] Please consider the following formatting changes --- PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx index 7a45fdfa525..ad3db4a7091 100644 --- a/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx +++ b/PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx @@ -220,8 +220,8 @@ struct AngularCorrelationsInJets { } if (doprocessRun2) { if (!(track.trackType() & o2::aod::track::Run2Track) //|| - //!(track.flags() & o2::aod::track::TPCrefit) || - //!(track.flags() & o2::aod::track::ITSrefit) + //!(track.flags() & o2::aod::track::TPCrefit) || + //!(track.flags() & o2::aod::track::ITSrefit) ) { return false; }