From 89fcaea19bf1b371dda8288fe8a38402865317ea Mon Sep 17 00:00:00 2001 From: Durgesh Bhatt Date: Sat, 30 May 2026 20:22:33 +0530 Subject: [PATCH] addedpt dependent cuts and configurables for qa plots --- PWGLF/Tasks/Resonances/CMakeLists.txt | 4 +- PWGLF/Tasks/Resonances/deltaAnalysis.cxx | 1116 ++++++++++++++++++++++ PWGLF/Tasks/Resonances/deltaanalysis.cxx | 613 ------------ 3 files changed, 1118 insertions(+), 615 deletions(-) create mode 100644 PWGLF/Tasks/Resonances/deltaAnalysis.cxx delete mode 100644 PWGLF/Tasks/Resonances/deltaanalysis.cxx diff --git a/PWGLF/Tasks/Resonances/CMakeLists.txt b/PWGLF/Tasks/Resonances/CMakeLists.txt index 1dd8a1015d8..74c48c0872a 100644 --- a/PWGLF/Tasks/Resonances/CMakeLists.txt +++ b/PWGLF/Tasks/Resonances/CMakeLists.txt @@ -94,8 +94,8 @@ o2physics_add_dpl_workflow(lambda1520spherocityanalysis PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(delta-analysis - SOURCES deltaanalysis.cxx +o2physics_add_dpl_workflow(deltaanalysis + SOURCES deltaAnalysis.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) diff --git a/PWGLF/Tasks/Resonances/deltaAnalysis.cxx b/PWGLF/Tasks/Resonances/deltaAnalysis.cxx new file mode 100644 index 00000000000..c2cc29ac62d --- /dev/null +++ b/PWGLF/Tasks/Resonances/deltaAnalysis.cxx @@ -0,0 +1,1116 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file deltaanalysis.cxx +/// \brief Delta(1232) resonance analysis via proton-pion invariant mass reconstruction with advance PID and background rejection cuts +/// \author Durgesh Bhatt + +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::constants::physics; + +namespace +{ +constexpr float massProton = o2::constants::physics::MassProton; +constexpr float massPion = o2::constants::physics::MassPionCharged; +} // namespace +namespace delta_analysis +{ +static constexpr int PdgProton{2212}; +static constexpr int PdgPion{211}; +static constexpr int PdgDeltaPlusPlus{2224}; +static constexpr int PdgDeltaZero{2114}; + +enum CentEstimator : int { + kFT0M = 0, + kFT0A = 1, + kFT0C = 2, + kFV0A = 3, + kNTPV = 4 +}; +} // namespace delta_analysis + +struct DeltaAnalysis { + + SliceCache cache; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + // FIX name/configurable: single space between type and name; member name == JSON key; lowerCamelCase + Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted |z-vertex| range [cm]"}; + Configurable applyOccupancyInTimeRangeCut{"applyOccupancyInTimeRangeCut", false, "Apply occupancy-in-time-range cut"}; + Configurable cfgOccupancyMin{"cfgOccupancyMin", 0, "Minimum track occupancy in time range"}; + Configurable cfgOccupancyMax{"cfgOccupancyMax", 9999, "Maximum track occupancy in time range"}; + + Configurable cfgCentralityEstimator{"cfgCentralityEstimator", 0, "Centrality estimator: 0=FT0M 1=FT0A 2=FT0C 3=FV0A 4=NTPV"}; + Configurable cfgCentMin{"cfgCentMin", 0.f, "Minimum centrality percentile"}; + Configurable cfgCentMax{"cfgCentMax", 100.f, "Maximum centrality percentile"}; + + Configurable cfgCutPt{"cfgCutPt", 0.2f, "Minimum pT of daughter track [GeV/c]"}; + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Maximum |eta| of daughter track"}; + // CHANGED: replaced single symmetric cfgCutY with independent cfgMinY / cfgMaxY + Configurable cfgMinY{"cfgMinY", -0.5f, "Minimum rapidity of reconstructed Delta"}; + Configurable cfgMaxY{"cfgMaxY", 0.5f, "Maximum rapidity of reconstructed Delta"}; + + Configurable cfgMinITSClusters{"cfgMinITSClusters", 5, "Minimum ITS clusters"}; + Configurable cfgMinTPCClusters{"cfgMinTPCClusters", 70, "Minimum TPC clusters found"}; + Configurable numberOfInvMassBins{"numberOfInvMassBins", 120, "Number of bins along the invariant mass axis"}; + Configurable cfgMinTPCCrossedRows{"cfgMinTPCCrossedRows", 70, "Minimum TPC crossed pad rows"}; + Configurable cfgMinCrossedRowsOverFindable{"cfgMinCrossedRowsOverFindable", 0.8f, "Minimum ratio crossedRows/findableClusters"}; + Configurable cfgMaxTPCSharedClusters{"cfgMaxTPCSharedClusters", 0, "Maximum TPC shared clusters"}; + Configurable cfgMaxTPCChi2NCl{"cfgMaxTPCChi2NCl", 4.0f, "Maximum TPC chi2/NCl"}; + Configurable cfgMaxITSChi2NCl{"cfgMaxITSChi2NCl", 36.0f, "Maximum ITS chi2/NCl"}; + Configurable requirePrimaryTrack{"requirePrimaryTrack", true, "Require isPrimaryTrack flag"}; + Configurable requireGlobalTrackNoDCA{"requireGlobalTrackNoDCA", true, "Require isGlobalTrackWoDCA flag"}; + Configurable requirePVContributor{"requirePVContributor", true, "Require PV-contributor flag"}; + + Configurable cfgCutDCAz{"cfgCutDCAz", 0.1f, "Maximum |DCAz| for all tracks [cm]"}; + Configurable> protonDCAPtBinEdges{"protonDCAPtBinEdges", {0.0f, 0.5f, 1.0f, 2.0f, 1000.f}, "Proton pT bin edges for DCAxy cut [GeV/c]"}; + Configurable> protonMaxDCAxyPerPtBin{"protonMaxDCAxyPerPtBin", {0.10f, 0.08f, 0.05f, 0.05f}, "Max |DCAxy| for proton per pT bin [cm]"}; + Configurable> pionDCAPtBinEdges{"pionDCAPtBinEdges", {0.0f, 0.5f, 1.0f, 2.0f, 1000.f}, "Pion pT bin edges for DCAxy cut [GeV/c]"}; + Configurable> pionMaxDCAxyPerPtBin{"pionMaxDCAxyPerPtBin", {0.20f, 0.15f, 0.10f, 0.08f}, "Max |DCAxy| for pion per pT bin [cm]"}; + + Configurable useTPCOnlyPID{"useTPCOnlyPID", false, "Use TPC-only PID (ignore TOF even if present)"}; + Configurable requireTOFForProton{"requireTOFForProton", false, "Require TOF signal for proton candidates"}; + Configurable requireTOFForPion{"requireTOFForPion", false, "Require TOF signal for pion candidates"}; + + Configurable minProtonMomentum{"minProtonMomentum", 0.f, "Minimum proton momentum for TOF PID [GeV/c]"}; + Configurable minTPCNSigmaProton{"minTPCNSigmaProton", -6.0f, "Minimum (lower bound) TPC nSigma for proton"}; + Configurable minTOFNSigmaProton{"minTOFNSigmaProton", -6.0f, "Minimum (lower bound) TOF nSigma for proton"}; + Configurable minCombinedNSigmaProton{"minCombinedNSigmaProton", -6.0f, "Minimum combined nSigma for proton (asymmetric mode)"}; + Configurable maxTPCNSigmaProton{"maxTPCNSigmaProton", 3.0, "Maximum |TPC nSigma| for proton"}; + Configurable combinedNSigmaCutProton{"combinedNSigmaCutProton", 3.0, "Circular TPC+TOF combined nSigma cut for proton. Negative = asymmetric mode."}; + Configurable> protonTPCPIDMomentumBins{"protonTPCPIDMomentumBins", {0.f, 0.5f, 0.7f, 0.8f}, "Proton TPC PID momentum bin edges [GeV/c]"}; + Configurable> protonTPCNSigmaCutPerBin{"protonTPCNSigmaCutPerBin", {5.f, 3.5f, 2.5f}, "Maximum TPC nSigma for proton per momentum bin"}; + Configurable> protonTOFPIDMomentumBins{"protonTOFPIDMomentumBins", {0.f, 999.f}, "Proton TOF PID momentum bin edges [GeV/c]"}; + Configurable> protonTOFNSigmaCutPerBin{"protonTOFNSigmaCutPerBin", {3.0f}, "Maximum TOF nSigma for proton per momentum bin"}; + + Configurable minPionMomentum{"minPionMomentum", 0.f, "Minimum pion momentum for TOF PID [GeV/c]"}; + Configurable minTPCNSigmaPion{"minTPCNSigmaPion", -6.0f, "Minimum (lower bound) TPC nSigma for pion"}; + Configurable minTOFNSigmaPion{"minTOFNSigmaPion", -6.0f, "Minimum (lower bound) TOF nSigma for pion"}; + Configurable minCombinedNSigmaPion{"minCombinedNSigmaPion", -6.0f, "Minimum combined nSigma for pion (asymmetric mode)"}; + Configurable maxTPCNSigmaPion{"maxTPCNSigmaPion", 3.0, "Maximum |TPC nSigma| for pion"}; + Configurable combinedNSigmaCutPion{"combinedNSigmaCutPion", 3.0, "Circular TPC+TOF combined nSigma cut for pion. Negative = asymmetric mode."}; + Configurable> pionTPCPIDMomentumBins{"pionTPCPIDMomentumBins", {0.f, 0.25f, 0.4f, 0.5f}, "Pion TPC PID momentum bin edges [GeV/c]"}; + Configurable> pionTPCNSigmaCutPerBin{"pionTPCNSigmaCutPerBin", {5.f, 3.5f, 2.5f}, "Maximum TPC nSigma for pion per momentum bin"}; + Configurable> pionTOFPIDMomentumBins{"pionTOFPIDMomentumBins", {0.f, 999.f}, "Pion TOF PID momentum bin edges [GeV/c]"}; + Configurable> pionTOFNSigmaCutPerBin{"pionTOFNSigmaCutPerBin", {3.0f}, "Maximum TOF nSigma for pion per momentum bin"}; + + Configurable tpcNSigmaVetoThreshold{"tpcNSigmaVetoThreshold", 3.0f, "Reject track if TPC nSigma of a competing species is below this value"}; + Configurable tofNSigmaVetoThreshold{"tofNSigmaVetoThreshold", 3.0f, "Reject track if TOF nSigma of a competing species is below this value"}; + + Configurable applyDeepAngleCut{"applyDeepAngleCut", false, "Apply minimum opening-angle cut (removes split-track background)"}; + Configurable deepAngleCutValue{"deepAngleCutValue", 0.04, "Minimum opening angle between proton and pion [rad]"}; + + Configurable cfgNoMixedEvents{"cfgNoMixedEvents", 5, "Number of mixed events per signal event"}; + + Configurable enableRotationalBackground{"enableRotationalBackground", false, "Compute rotational background by rotating pion phi near PI"}; + Configurable numberOfRotations{"numberOfRotations", 10, "Number of pion-phi rotations for background"}; + Configurable rotationAngleWindow{"rotationAngleWindow", 6.f, "Pion rotated by angles within PI +/- PI/rotationAngleWindow"}; + + ConfigurableAxis cfgPtAxis{"cfgPtAxis", {VARIABLE_WIDTH, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4.0, 5.0, 7.0, 10.0}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis cfgCentAxis{"cfgCentAxis", {VARIABLE_WIDTH, 0.f, 10.f, 20.f, 30.f, 40.f, 50.f, 60.f, 70.f, 80.f, 90.f, 100.f}, "Centrality (%)"}; + ConfigurableAxis cfgVtxAxis{"cfgVtxAxis", {VARIABLE_WIDTH, -12.f, -10.f, -9.f, -8.f, -7.f, -6.f, -5.f, -4.f, -3.f, -2.f, -1.f, 0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 12.f}, "Vertex z [cm]"}; + ConfigurableAxis cfgRapAxis{"cfgRapAxis", {20, -1.0, 1.0}, "Rapidity y"}; + + std::vector mProtonTPCMomBins{}; + std::vector mProtonTPCNSigCuts{}; + std::vector mProtonTOFMomBins{}; + std::vector mProtonTOFNSigCuts{}; + std::vector mPionTPCMomBins{}; + std::vector mPionTPCNSigCuts{}; + std::vector mPionTOFMomBins{}; + std::vector mPionTOFNSigCuts{}; + std::vector mProtonDCAPtEdges{}; + std::vector mProtonMaxDCAxy{}; + std::vector mPionDCAPtEdges{}; + std::vector mPionMaxDCAxy{}; + + void init(InitContext const&) + { + auto chkSz = [](auto const& cuts, auto const& bins, const char* label) { + const auto& cutsVec = cuts.value; + const auto& binsVec = bins.value; + if (cutsVec.size() != binsVec.size() - 1) { + LOG(fatal) << "PID/DCA bin-size mismatch for: " << label + << " (cuts=" << cutsVec.size() + << ", bins=" << binsVec.size() << ")"; + } + }; + chkSz(protonTPCNSigmaCutPerBin, protonTPCPIDMomentumBins, "proton TPC PID"); + chkSz(protonTOFNSigmaCutPerBin, protonTOFPIDMomentumBins, "proton TOF PID"); + chkSz(pionTPCNSigmaCutPerBin, pionTPCPIDMomentumBins, "pion TPC PID"); + chkSz(pionTOFNSigmaCutPerBin, pionTOFPIDMomentumBins, "pion TOF PID"); + chkSz(protonMaxDCAxyPerPtBin, protonDCAPtBinEdges, "proton DCA"); + chkSz(pionMaxDCAxyPerPtBin, pionDCAPtBinEdges, "pion DCA"); + + mProtonTPCMomBins = protonTPCPIDMomentumBins.value; + mProtonTPCNSigCuts = protonTPCNSigmaCutPerBin.value; + mProtonTOFMomBins = protonTOFPIDMomentumBins.value; + mProtonTOFNSigCuts = protonTOFNSigmaCutPerBin.value; + mPionTPCMomBins = pionTPCPIDMomentumBins.value; + mPionTPCNSigCuts = pionTPCNSigmaCutPerBin.value; + mPionTOFMomBins = pionTOFPIDMomentumBins.value; + mPionTOFNSigCuts = pionTOFNSigmaCutPerBin.value; + mProtonDCAPtEdges = protonDCAPtBinEdges.value; + mProtonMaxDCAxy = protonMaxDCAxyPerPtBin.value; + mPionDCAPtEdges = pionDCAPtBinEdges.value; + mPionMaxDCAxy = pionMaxDCAxyPerPtBin.value; + + const AxisSpec ptAxis{200, 0., 10., "p_{T} (GeV/c)"}; + const AxisSpec massAxis{numberOfInvMassBins, 0.8, 1.8, "M_{inv} (GeV/#it{c}^{2})"}; + const AxisSpec centAxis{cfgCentAxis, "Centrality (%)"}; + const AxisSpec vtxAxis{cfgVtxAxis, "Vertex z [cm]"}; + const AxisSpec rapAxis{cfgRapAxis, "Rapidity y"}; + const AxisSpec nSigmaTPCaxis{100, -10., 10., "n#sigma^{TPC}"}; + const AxisSpec nSigmaTOFaxis{100, -10., 10., "n#sigma^{TOF}"}; + const AxisSpec momentumAxis{200, 0., 10., "p (GeV/#it{c})"}; + const AxisSpec ptForPIDAxis{200, 0., 10., "p_{T} (GeV/#it{c})"}; + const AxisSpec dcaXYaxis{200, -0.2f, 0.2f, "DCA_{xy} (cm)"}; + const AxisSpec dcaZaxis{200, -1.2f, 1.2f, "DCA_{z} (cm)"}; + const AxisSpec tpcClusAxis{200, 0, 200, "N_{clusters}^{TPC}"}; + const AxisSpec tpcRowsAxis{200, 0, 200, "Crossed rows^{TPC}"}; + const AxisSpec occupancyAxis{200, 0, 10000, "Track occupancy [-40, 100] ns"}; + const AxisSpec etaAxis{100, -1.0, 1.0, "#eta"}; + const AxisSpec openAngleAxis{100, 0., o2::constants::math::PI, "Opening angle (rad)"}; + + histos.add("Event/hVtxZ", "Vertex z; z (cm)", kTH1F, {{400, -20., 20.}}); + histos.add("Event/hNcontributor", "PV contributors; N", kTH1F, {{2001, -0.5f, 2000.5f}}); + histos.add("Event/hCentrality", "Centrality", kTH1F, {centAxis}); + histos.add("Event/hOccupancy", "Occupancy in time range", kTH1F, {occupancyAxis}); + + histos.add("CentQA/hCentralityVsVtxZ", "Centrality vs vertex z", kTH2F, {vtxAxis, centAxis}); + histos.add("CentQA/hCentralityVsOccupancy", "Centrality vs occupancy", kTH2F, {occupancyAxis, centAxis}); + histos.add("CentQA/hEventCountVsCentrality", "Event count vs centrality", kTH1F, {centAxis}); + + histos.add("OccupancyQA/hOccupancyVsCentralityBefore", "Occupancy vs centrality (before cut)", kTH2F, {centAxis, occupancyAxis}); + histos.add("OccupancyQA/hOccupancyVsVtxZBefore", "Occupancy vs vertex z (before cut)", kTH2F, {vtxAxis, occupancyAxis}); + histos.add("OccupancyQA/hOccupancyVsCentralityAfter", "Occupancy vs centrality (after cut)", kTH2F, {centAxis, occupancyAxis}); + histos.add("OccupancyQA/hOccupancyVsVtxZAfter", "Occupancy vs vertex z (after cut)", kTH2F, {vtxAxis, occupancyAxis}); + + histos.add("QAbefore/Proton/tpcNSigmaVsMomentum", "TPC n#sigma proton vs p (before cuts)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("QAbefore/Proton/tofNSigmaVsMomentum", "TOF n#sigma proton vs p (before cuts)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("QAbefore/Proton/tofNSigmaVsTPCNSigma", "TOF vs TPC n#sigma proton (before cuts)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + histos.add("QAbefore/Pion/tpcNSigmaVsMomentum", "TPC n#sigma pion vs p (before cuts)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("QAbefore/Pion/tofNSigmaVsMomentum", "TOF n#sigma pion vs p (before cuts)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("QAbefore/Pion/tofNSigmaVsTPCNSigma", "TOF vs TPC n#sigma pion (before cuts)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + + histos.add("QAafter/Proton/dcaXYvsPt", "Proton DCA_{xy} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaXYaxis}); + histos.add("QAafter/Proton/dcaZvsPt", "Proton DCA_{z} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaZaxis}); + histos.add("QAafter/Proton/tpcNSigmaVsMomentum", "Proton TPC n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("QAafter/Proton/tpcNSigmaVsPt", "Proton TPC n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("QAafter/Proton/tpcNSigmaVsCentrality", "Proton TPC n#sigma vs centrality", kTH2F, {centAxis, nSigmaTPCaxis}); + histos.add("QAafter/Proton/tofNSigmaVsMomentum", "Proton TOF n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("QAafter/Proton/tofNSigmaVsPt", "Proton TOF n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTOFaxis}); + histos.add("QAafter/Proton/tofNSigmaVsCentrality", "Proton TOF n#sigma vs centrality", kTH2F, {centAxis, nSigmaTOFaxis}); + histos.add("QAafter/Proton/tofNSigmaVsTPCNSigma", "Proton TOF vs TPC n#sigma (after cuts)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + histos.add("QAafter/Proton/tpcNSigmaPionContamVsPt", "Proton track: TPC n#sigma pion contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("QAafter/Proton/tpcNSigmaKaonContamVsPt", "Proton track: TPC n#sigma kaon contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("QAafter/Proton/tofNSigmaPionContamVsMomentum", "Proton track: TOF n#sigma pion contamination", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("QAafter/Proton/tpcCrossedRowsVsPt", "Proton TPC crossed rows vs p_{T}", kTH2F, {ptForPIDAxis, tpcRowsAxis}); + histos.add("QAafter/Proton/tpcClustersFoundVsPt", "Proton TPC clusters found vs p_{T}", kTH2F, {ptForPIDAxis, tpcClusAxis}); + histos.add("QAafter/Proton/dcaXYdist", "Proton DCA_{xy} distribution (fine bins)", kTH1F, {dcaXYaxis}); + histos.add("QAafter/Proton/dcaZdist", "Proton DCA_{z} distribution (fine bins)", kTH1F, {dcaZaxis}); + + histos.add("QAafter/Pion/dcaXYvsPt", "Pion DCA_{xy} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaXYaxis}); + histos.add("QAafter/Pion/dcaZvsPt", "Pion DCA_{z} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaZaxis}); + histos.add("QAafter/Pion/tpcNSigmaVsMomentum", "Pion TPC n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("QAafter/Pion/tpcNSigmaVsPt", "Pion TPC n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("QAafter/Pion/tpcNSigmaVsCentrality", "Pion TPC n#sigma vs centrality", kTH2F, {centAxis, nSigmaTPCaxis}); + histos.add("QAafter/Pion/tofNSigmaVsMomentum", "Pion TOF n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("QAafter/Pion/tofNSigmaVsPt", "Pion TOF n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTOFaxis}); + histos.add("QAafter/Pion/tofNSigmaVsCentrality", "Pion TOF n#sigma vs centrality", kTH2F, {centAxis, nSigmaTOFaxis}); + histos.add("QAafter/Pion/tofNSigmaVsTPCNSigma", "Pion TOF vs TPC n#sigma (after cuts)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + histos.add("QAafter/Pion/tpcNSigmaProtonContamVsPt", "Pion track: TPC n#sigma proton contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("QAafter/Pion/tpcNSigmaKaonContamVsPt", "Pion track: TPC n#sigma kaon contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("QAafter/Pion/tofNSigmaProtonContamVsMomentum", "Pion track: TOF n#sigma proton contamination", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("QAafter/Pion/tpcCrossedRowsVsPt", "Pion TPC crossed rows vs p_{T}", kTH2F, {ptForPIDAxis, tpcRowsAxis}); + histos.add("QAafter/Pion/tpcClustersFoundVsPt", "Pion TPC clusters found vs p_{T}", kTH2F, {ptForPIDAxis, tpcClusAxis}); + histos.add("QAafter/Pion/dcaXYdist", "Pion DCA_{xy} distribution (fine bins)", kTH1F, {dcaXYaxis}); + histos.add("QAafter/Pion/dcaZdist", "Pion DCA_{z} distribution (fine bins)", kTH1F, {dcaZaxis}); + + histos.add("QAChecks/Pair/hOpeningAngleBefore", "Opening angle (before cuts)", kTH1F, {openAngleAxis}); + histos.add("QAChecks/Pair/hOpeningAngleAfter", "Opening angle (after cuts)", kTH1F, {openAngleAxis}); + + histos.add("Analysis/hDeltaPlusPlusInvMass", "#Delta^{++} (#rightarrow p + #pi^{+}) invariant mass", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaPlusPlusInvMass", "#bar{#Delta}^{++} (#rightarrow #bar{p} + #pi^{-}) invariant mass", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hDeltaZeroInvMass", "#Delta^{0} (#rightarrow p + #pi^{-}) invariant mass", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaZeroInvMass", "#bar{#Delta}^{0} (#rightarrow #bar{p} + #pi^{+}) invariant mass", kTH2F, {ptAxis, massAxis}); + + if (enableRotationalBackground) { + histos.add("Analysis/hDeltaPlusPlusInvMassRot", "#Delta^{++} invariant mass - rotational background", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaPlusPlusInvMassRot", "#bar{#Delta}^{++} invariant mass - rotational background", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hDeltaZeroInvMassRot", "#Delta^{0} invariant mass - rotational background", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaZeroInvMassRot", "#bar{#Delta}^{0} invariant mass - rotational background", kTH2F, {ptAxis, massAxis}); + } + + if (doprocessMixedEvent) { + histos.add("Analysis/hDeltaPlusPlusInvMassEM", "#Delta^{++} invariant mass - event mixing", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaPlusPlusInvMassEM", "#bar{#Delta}^{++} invariant mass - event mixing", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hDeltaZeroInvMassEM", "#Delta^{0} invariant mass - event mixing", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaZeroInvMassEM", "#bar{#Delta}^{0} invariant mass - event mixing", kTH2F, {ptAxis, massAxis}); + } + + histos.add("THnSparse/hDeltaPlusPlus", "THnSparse #Delta^{++} same-event", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaPlusPlus", "THnSparse #bar{#Delta}^{++} same-event", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hDeltaZero", "THnSparse #Delta^{0} same-event", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaZero", "THnSparse #bar{#Delta}^{0} same-event", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + + if (enableRotationalBackground) { + histos.add("THnSparse/hDeltaPlusPlusRot", "THnSparse #Delta^{++} rotational background", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaPlusPlusRot", "THnSparse #bar{#Delta}^{++} rotational background", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hDeltaZeroRot", "THnSparse #Delta^{0} rotational background", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaZeroRot", "THnSparse #bar{#Delta}^{0} rotational background", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + } + + if (doprocessMixedEvent) { + histos.add("THnSparse/hDeltaPlusPlusEM", "THnSparse #Delta^{++} event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaPlusPlusEM", "THnSparse #bar{#Delta}^{++} event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hDeltaZeroEM", "THnSparse #Delta^{0} event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaZeroEM", "THnSparse #bar{#Delta}^{0} event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + } + + if (doprocessMC) { + histos.add("THnSparse/hDeltaPlusPlusMC", "THnSparse #Delta^{++} MC reconstructed", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaPlusPlusMC", "THnSparse #bar{#Delta}^{++} MC reconstructed", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hDeltaZeroMC", "THnSparse #Delta^{0} MC reconstructed", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaZeroMC", "THnSparse #bar{#Delta}^{0} MC reconstructed", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hDeltaPlusPlusGen", "THnSparse #Delta^{++} generated", kTHnSparseF, {ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaPlusPlusGen", "THnSparse #bar{#Delta}^{++} generated", kTHnSparseF, {ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hDeltaZeroGen", "THnSparse #Delta^{0} generated", kTHnSparseF, {ptAxis, centAxis, rapAxis}); + histos.add("THnSparse/hAntiDeltaZeroGen", "THnSparse #bar{#Delta}^{0} generated", kTHnSparseF, {ptAxis, centAxis, rapAxis}); + + histos.add("Analysis/hDeltaPlusPlusInvMassMC", "#Delta^{++} invariant mass (MC truth-matched)", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaPlusPlusInvMassMC", "#bar{#Delta}^{++} invariant mass (MC truth-matched)", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hDeltaZeroInvMassMC", "#Delta^{0} invariant mass (MC truth-matched)", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaZeroInvMassMC", "#bar{#Delta}^{0} invariant mass (MC truth-matched)", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hDeltaPlusPlusInvMassGen", "#Delta^{++} invariant mass - generated", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaPlusPlusInvMassGen", "#bar{#Delta}^{++} invariant mass - generated", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hDeltaZeroInvMassGen", "#Delta^{0} invariant mass - generated", kTH2F, {ptAxis, massAxis}); + histos.add("Analysis/hAntiDeltaZeroInvMassGen", "#bar{#Delta}^{0} invariant mass - generated", kTH2F, {ptAxis, massAxis}); + + histos.add("QAChecks/hRecProtonDeltaPlusPlus", "Rec proton from #Delta^{++}", kTH1F, {ptAxis}); + histos.add("QAChecks/hRecProtonAntiDeltaPlusPlus", "Rec proton from #bar{#Delta}^{++}", kTH1F, {ptAxis}); + histos.add("QAChecks/hRecProtonDeltaZero", "Rec proton from #Delta^{0}", kTH1F, {ptAxis}); + histos.add("QAChecks/hRecProtonAntiDeltaZero", "Rec proton from #bar{#Delta}^{0}", kTH1F, {ptAxis}); + histos.add("QAChecks/hRecPionDeltaPlusPlus", "Rec pion from #Delta^{++}", kTH1F, {ptAxis}); + histos.add("QAChecks/hRecPionAntiDeltaPlusPlus", "Rec pion from #bar{#Delta}^{++}", kTH1F, {ptAxis}); + histos.add("QAChecks/hRecPionDeltaZero", "Rec pion from #Delta^{0}", kTH1F, {ptAxis}); + histos.add("QAChecks/hRecPionAntiDeltaZero", "Rec pion from #bar{#Delta}^{0}", kTH1F, {ptAxis}); + histos.add("QAChecks/hGenProtonDeltaPlusPlus", "Gen proton from #Delta^{++}", kTH1F, {ptAxis}); + histos.add("QAChecks/hGenProtonAntiDeltaPlusPlus", "Gen proton from #bar{#Delta}^{++}", kTH1F, {ptAxis}); + histos.add("QAChecks/hGenProtonDeltaZero", "Gen proton from #Delta^{0}", kTH1F, {ptAxis}); + histos.add("QAChecks/hGenProtonAntiDeltaZero", "Gen proton from #bar{#Delta}^{0}", kTH1F, {ptAxis}); + histos.add("QAChecks/hGenPionDeltaPlusPlus", "Gen pion from #Delta^{++}", kTH1F, {ptAxis}); + histos.add("QAChecks/hGenPionAntiDeltaPlusPlus", "Gen pion from #bar{#Delta}^{++}", kTH1F, {ptAxis}); + histos.add("QAChecks/hGenPionDeltaZero", "Gen pion from #Delta^{0}", kTH1F, {ptAxis}); + histos.add("QAChecks/hGenPionAntiDeltaZero", "Gen pion from #bar{#Delta}^{0}", kTH1F, {ptAxis}); + } + } // end init() + + template + float getCentrality(CollisionType const& collision) + { + switch (cfgCentralityEstimator) { + case delta_analysis::kFT0M: + return collision.centFT0M(); + case delta_analysis::kFT0A: + return collision.centFT0A(); + case delta_analysis::kFT0C: + return collision.centFT0C(); + case delta_analysis::kFV0A: + return collision.centFV0A(); + case delta_analysis::kNTPV: + return collision.centNTPV(); + default: + return collision.centFT0M(); + } + } + + template + bool passesEventSelection(CollisionType const& collision) + { + if (!collision.sel8()) + return false; + if (std::abs(collision.posZ()) > cfgCutVertex) + return false; + if (applyOccupancyInTimeRangeCut) { + const int occ = collision.trackOccupancyInTimeRange(); + if (occ < cfgOccupancyMin || occ > cfgOccupancyMax) + return false; + } + const float cent = getCentrality(collision); + if (cent < cfgCentMin || cent > cfgCentMax) + return false; + return true; + } + + template + bool passesBasicTrackSelection(TrackType const& track) + { + if (track.itsNCls() < cfgMinITSClusters) + return false; + if (track.tpcNClsShared() > cfgMaxTPCSharedClusters) + return false; + if (track.tpcNClsFound() < cfgMinTPCClusters) + return false; + if (track.tpcNClsCrossedRows() < cfgMinTPCCrossedRows) + return false; + if (track.tpcNClsCrossedRows() < cfgMinCrossedRowsOverFindable * track.tpcNClsFindable()) + return false; + if (track.tpcChi2NCl() > cfgMaxTPCChi2NCl) + return false; + if (track.itsChi2NCl() > cfgMaxITSChi2NCl) + return false; + if (requirePrimaryTrack && !track.isPrimaryTrack()) + return false; + if (requireGlobalTrackNoDCA && !track.isGlobalTrackWoDCA()) + return false; + if (requirePVContributor && !track.isPVContributor()) + return false; + return true; + } + + template + bool passesProtonDCASelection(TrackType const& track) + { + const int nBins = static_cast(mProtonDCAPtEdges.size()) - 1; + const float pt = track.pt(); + bool passed = false; + for (int i = 0; i < nBins; ++i) { + if (pt >= mProtonDCAPtEdges[i] && pt < mProtonDCAPtEdges[i + 1] && + std::abs(track.dcaXY()) < mProtonMaxDCAxy[i]) { + passed = true; + break; + } + } + return passed && (std::abs(track.dcaZ()) < cfgCutDCAz); + } + + template + bool passesPionDCASelection(TrackType const& track) + { + const int nBins = static_cast(mPionDCAPtEdges.size()) - 1; + const float pt = track.pt(); + bool passed = false; + for (int i = 0; i < nBins; ++i) { + if (pt >= mPionDCAPtEdges[i] && pt < mPionDCAPtEdges[i + 1] && + std::abs(track.dcaXY()) < mPionMaxDCAxy[i]) { + passed = true; + break; + } + } + return passed && (std::abs(track.dcaZ()) < cfgCutDCAz); + } + + template + bool passesProtonPID(TrackType const& track, float totalMomentum) + { + bool tpcPassed{false}, tofPassed{false}; + const int nTPCBins = static_cast(mProtonTPCMomBins.size()); + const int nTOFBins = static_cast(mProtonTOFMomBins.size()); + const float tpcNSigPi = std::abs(track.tpcNSigmaPi()); + const float tpcNSigPr = std::abs(track.tpcNSigmaPr()); + const float tofNSigPi = std::abs(track.tofNSigmaPi()); + const float tofNSigPr = std::abs(track.tofNSigmaPr()); + const float combinedNSigPr = tpcNSigPr * tpcNSigPr + tofNSigPr * tofNSigPr; + const float combinedNSigPi = tpcNSigPi * tpcNSigPi + tofNSigPi * tofNSigPi; + const float circularCutSq = static_cast(combinedNSigmaCutProton * combinedNSigmaCutProton); + + const float circularVetoCutSq = tpcNSigmaVetoThreshold * tpcNSigmaVetoThreshold + tofNSigmaVetoThreshold * tofNSigmaVetoThreshold; + + if (!useTPCOnlyPID && track.hasTOF()) { + if (combinedNSigmaCutProton < 0 && totalMomentum >= minProtonMomentum) { + if (track.tofNSigmaPr() < minTOFNSigmaProton) + return false; + for (int i = 0; i < nTOFBins - 1; ++i) { + if (totalMomentum >= mProtonTOFMomBins[i] && totalMomentum < mProtonTOFMomBins[i + 1] && + tofNSigPr < mProtonTOFNSigCuts[i] && tofNSigPi > tofNSigmaVetoThreshold) { + tofPassed = true; + break; + } + } + if (track.tpcNSigmaPr() < minCombinedNSigmaProton) + return false; + if (tpcNSigPr < static_cast(maxTPCNSigmaProton) && tpcNSigPi > tpcNSigmaVetoThreshold) + tpcPassed = true; + } else if (combinedNSigmaCutProton > 0 && totalMomentum >= minProtonMomentum) { + if (combinedNSigPr < circularCutSq && combinedNSigPi > circularVetoCutSq) { + tpcPassed = true; + tofPassed = true; + } + } + if (totalMomentum < minProtonMomentum && tpcNSigPr < static_cast(maxTPCNSigmaProton)) { + tpcPassed = true; + tofPassed = true; + } + } else { + tofPassed = true; + if (track.tpcNSigmaPr() < minTPCNSigmaProton) + return false; + for (int i = 0; i < nTPCBins - 1; ++i) { + if (totalMomentum >= mProtonTPCMomBins[i] && totalMomentum < mProtonTPCMomBins[i + 1] && + tpcNSigPr < mProtonTPCNSigCuts[i] && tpcNSigPi > tpcNSigmaVetoThreshold) { + tpcPassed = true; + break; + } + } + } + return tpcPassed && tofPassed; + } + + template + bool passesPionPID(TrackType const& track, float totalMomentum) + { + bool tpcPassed{false}, tofPassed{false}; + const int nTPCBins = static_cast(mPionTPCMomBins.size()); + const int nTOFBins = static_cast(mPionTOFMomBins.size()); + const float tpcNSigPi = std::abs(track.tpcNSigmaPi()); + const float tpcNSigPr = std::abs(track.tpcNSigmaPr()); + const float tofNSigPi = std::abs(track.tofNSigmaPi()); + const float tofNSigPr = std::abs(track.tofNSigmaPr()); + const float combinedNSigPi = tpcNSigPi * tpcNSigPi + tofNSigPi * tofNSigPi; + const float combinedNSigPr = tpcNSigPr * tpcNSigPr + tofNSigPr * tofNSigPr; + const float circularCutSq = static_cast(combinedNSigmaCutPion * combinedNSigmaCutPion); + const float circularVetoCutSq = tpcNSigmaVetoThreshold * tpcNSigmaVetoThreshold + tofNSigmaVetoThreshold * tofNSigmaVetoThreshold; + + if (!useTPCOnlyPID && track.hasTOF()) { + if (combinedNSigmaCutPion < 0 && totalMomentum >= minPionMomentum) { + if (track.tofNSigmaPi() < minTOFNSigmaPion) + return false; + for (int i = 0; i < nTOFBins - 1; ++i) { + if (totalMomentum >= mPionTOFMomBins[i] && totalMomentum < mPionTOFMomBins[i + 1] && + tofNSigPi < mPionTOFNSigCuts[i] && tofNSigPr > tofNSigmaVetoThreshold) { + tofPassed = true; + break; + } + } + if (track.tpcNSigmaPi() < minCombinedNSigmaPion) + return false; + if (tpcNSigPi < static_cast(maxTPCNSigmaPion) && tpcNSigPr > tpcNSigmaVetoThreshold) + tpcPassed = true; + } else if (combinedNSigmaCutPion > 0 && totalMomentum >= minPionMomentum) { + if (combinedNSigPi < circularCutSq && combinedNSigPr > circularVetoCutSq) { + tpcPassed = true; + tofPassed = true; + } + } + if (totalMomentum < minPionMomentum && tpcNSigPi < static_cast(maxTPCNSigmaPion)) { + tpcPassed = true; + tofPassed = true; + } + } else { + tofPassed = true; + if (track.tpcNSigmaPi() < minTPCNSigmaPion) + return false; + for (int i = 0; i < nTPCBins - 1; ++i) { + if (totalMomentum >= mPionTPCMomBins[i] && totalMomentum < mPionTPCMomBins[i + 1] && + tpcNSigPi < mPionTPCNSigCuts[i] && tpcNSigPr > tpcNSigmaVetoThreshold) { + tpcPassed = true; + break; + } + } + } + return tpcPassed && tofPassed; + } + + template + void fillQAProton(TrackType const& track, float totalMomentum, float centralityPercent) + { + const float pt = track.pt(); + const float tpcNSigPr = track.tpcNSigmaPr(); + histos.fill(HIST("QAafter/Proton/dcaXYvsPt"), pt, track.dcaXY()); + histos.fill(HIST("QAafter/Proton/dcaZvsPt"), pt, track.dcaZ()); + histos.fill(HIST("QAafter/Proton/dcaXYdist"), track.dcaXY()); + histos.fill(HIST("QAafter/Proton/dcaZdist"), track.dcaZ()); + histos.fill(HIST("QAafter/Proton/tpcNSigmaVsMomentum"), totalMomentum, tpcNSigPr); + histos.fill(HIST("QAafter/Proton/tpcNSigmaVsPt"), pt, tpcNSigPr); + histos.fill(HIST("QAafter/Proton/tpcNSigmaVsCentrality"), centralityPercent, tpcNSigPr); + histos.fill(HIST("QAafter/Proton/tpcNSigmaPionContamVsPt"), pt, track.tpcNSigmaPi()); + histos.fill(HIST("QAafter/Proton/tpcNSigmaKaonContamVsPt"), pt, track.tpcNSigmaKa()); + histos.fill(HIST("QAafter/Proton/tpcCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); + histos.fill(HIST("QAafter/Proton/tpcClustersFoundVsPt"), pt, track.tpcNClsFound()); + if (!useTPCOnlyPID && track.hasTOF()) { + const float tofNSigPr = track.tofNSigmaPr(); + histos.fill(HIST("QAafter/Proton/tofNSigmaVsMomentum"), totalMomentum, tofNSigPr); + histos.fill(HIST("QAafter/Proton/tofNSigmaVsPt"), pt, tofNSigPr); + histos.fill(HIST("QAafter/Proton/tofNSigmaVsCentrality"), centralityPercent, tofNSigPr); + histos.fill(HIST("QAafter/Proton/tofNSigmaVsTPCNSigma"), tpcNSigPr, tofNSigPr); + histos.fill(HIST("QAafter/Proton/tofNSigmaPionContamVsMomentum"), totalMomentum, track.tofNSigmaPi()); + } + } + + template + void fillQAPion(TrackType const& track, float totalMomentum, float centralityPercent) + { + const float pt = track.pt(); + const float tpcNSigPi = track.tpcNSigmaPi(); + histos.fill(HIST("QAafter/Pion/dcaXYvsPt"), pt, track.dcaXY()); + histos.fill(HIST("QAafter/Pion/dcaZvsPt"), pt, track.dcaZ()); + histos.fill(HIST("QAafter/Pion/dcaXYdist"), track.dcaXY()); + histos.fill(HIST("QAafter/Pion/dcaZdist"), track.dcaZ()); + histos.fill(HIST("QAafter/Pion/tpcNSigmaVsMomentum"), totalMomentum, tpcNSigPi); + histos.fill(HIST("QAafter/Pion/tpcNSigmaVsPt"), pt, tpcNSigPi); + histos.fill(HIST("QAafter/Pion/tpcNSigmaVsCentrality"), centralityPercent, tpcNSigPi); + histos.fill(HIST("QAafter/Pion/tpcNSigmaProtonContamVsPt"), pt, track.tpcNSigmaPr()); + histos.fill(HIST("QAafter/Pion/tpcNSigmaKaonContamVsPt"), pt, track.tpcNSigmaKa()); + histos.fill(HIST("QAafter/Pion/tpcCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); + histos.fill(HIST("QAafter/Pion/tpcClustersFoundVsPt"), pt, track.tpcNClsFound()); + if (!useTPCOnlyPID && track.hasTOF()) { + const float tofNSigPi = track.tofNSigmaPi(); + histos.fill(HIST("QAafter/Pion/tofNSigmaVsMomentum"), totalMomentum, tofNSigPi); + histos.fill(HIST("QAafter/Pion/tofNSigmaVsPt"), pt, tofNSigPi); + histos.fill(HIST("QAafter/Pion/tofNSigmaVsCentrality"), centralityPercent, tofNSigPi); + histos.fill(HIST("QAafter/Pion/tofNSigmaVsTPCNSigma"), tpcNSigPi, tofNSigPi); + histos.fill(HIST("QAafter/Pion/tofNSigmaProtonContamVsMomentum"), totalMomentum, track.tofNSigmaPr()); + } + } + + void fillDeltaHistogramSameEvent(int protonSign, int pionSign, float pairPt, float pairMass, float centrality, float rapidity) + { + if (protonSign > 0) { + if (pionSign > 0) { + histos.fill(HIST("Analysis/hDeltaPlusPlusInvMass"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hDeltaPlusPlus"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("Analysis/hDeltaZeroInvMass"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hDeltaZero"), pairMass, pairPt, centrality, rapidity); + } + } else { + if (pionSign < 0) { + histos.fill(HIST("Analysis/hAntiDeltaPlusPlusInvMass"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hAntiDeltaPlusPlus"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("Analysis/hAntiDeltaZeroInvMass"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hAntiDeltaZero"), pairMass, pairPt, centrality, rapidity); + } + } + } + + void fillDeltaHistogramMixedEvent(int protonSign, int pionSign, float pairPt, float pairMass, float centrality, float rapidity) + { + if (protonSign > 0) { + if (pionSign > 0) { + histos.fill(HIST("Analysis/hDeltaPlusPlusInvMassEM"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hDeltaPlusPlusEM"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("Analysis/hDeltaZeroInvMassEM"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hDeltaZeroEM"), pairMass, pairPt, centrality, rapidity); + } + } else { + if (pionSign < 0) { + histos.fill(HIST("Analysis/hAntiDeltaPlusPlusInvMassEM"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hAntiDeltaPlusPlusEM"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("Analysis/hAntiDeltaZeroInvMassEM"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hAntiDeltaZeroEM"), pairMass, pairPt, centrality, rapidity); + } + } + } + + void fillDeltaHistogramMC(int protonSign, int pionSign, float pairPt, float pairMass, float centrality, float rapidity) + { + if (protonSign > 0) { + if (pionSign > 0) { + histos.fill(HIST("Analysis/hDeltaPlusPlusInvMassMC"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hDeltaPlusPlusMC"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("Analysis/hDeltaZeroInvMassMC"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hDeltaZeroMC"), pairMass, pairPt, centrality, rapidity); + } + } else { + if (pionSign < 0) { + histos.fill(HIST("Analysis/hAntiDeltaPlusPlusInvMassMC"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hAntiDeltaPlusPlusMC"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("Analysis/hAntiDeltaZeroInvMassMC"), pairPt, pairMass); + histos.fill(HIST("THnSparse/hAntiDeltaZeroMC"), pairMass, pairPt, centrality, rapidity); + } + } + } + + void fillRotationalBackground(int protonSign, int pionSign, float pxProton, float pyProton, float pzProton, float pionPhi, float pionPt, float pzPion, float centrality) + { + if (numberOfRotations <= 0) + return; + const float weight = 1.f / static_cast(numberOfRotations); + const float rotWindowHalf = o2::constants::math::PI / static_cast(rotationAngleWindow); + + for (int iRot = 0; iRot < numberOfRotations; ++iRot) { + float rotAngle; + if (numberOfRotations == 1) { + rotAngle = o2::constants::math::PI; + } else { + rotAngle = (o2::constants::math::PI - rotWindowHalf) + + static_cast(iRot) * (2.f * rotWindowHalf / static_cast(numberOfRotations - 1)); + } + const float newPhi = RecoDecay::constrainAngle(pionPhi + rotAngle, 0.f); + const float pxPionRot = pionPt * std::cos(newPhi); + const float pyPionRot = pionPt * std::sin(newPhi); + const std::array, 2> rotMomenta = { + std::array{pxProton, pyProton, pzProton}, + std::array{pxPionRot, pyPionRot, pzPion}}; + const float rotMass = RecoDecay::m(rotMomenta, std::array{massProton, massPion}); + const float rotPt = RecoDecay::pt(std::array{pxProton + pxPionRot, pyProton + pyPionRot}); + const float rotY = RecoDecay::y( + std::array{pxProton + pxPionRot, pyProton + pyPionRot, pzProton + pzPion}, rotMass); + if (rotY < cfgMinY || rotY > cfgMaxY) + continue; + + if (protonSign > 0) { + if (pionSign > 0) { + histos.fill(HIST("Analysis/hDeltaPlusPlusInvMassRot"), rotPt, rotMass, weight); + histos.fill(HIST("THnSparse/hDeltaPlusPlusRot"), rotMass, rotPt, centrality, rotY, weight); + } else { + histos.fill(HIST("Analysis/hDeltaZeroInvMassRot"), rotPt, rotMass, weight); + histos.fill(HIST("THnSparse/hDeltaZeroRot"), rotMass, rotPt, centrality, rotY, weight); + } + } else { + if (pionSign < 0) { + histos.fill(HIST("Analysis/hAntiDeltaPlusPlusInvMassRot"), rotPt, rotMass, weight); + histos.fill(HIST("THnSparse/hAntiDeltaPlusPlusRot"), rotMass, rotPt, centrality, rotY, weight); + } else { + histos.fill(HIST("Analysis/hAntiDeltaZeroInvMassRot"), rotPt, rotMass, weight); + histos.fill(HIST("THnSparse/hAntiDeltaZeroRot"), rotMass, rotPt, centrality, rotY, weight); + } + } + } + } + + void fillPairQABefore(float pxPr, float pyPr, float pzPr, float pxPi, float pyPi, float pzPi, float momPr, float momPi) + { + const float cosAngle = std::clamp((pxPr * pxPi + pyPr * pyPi + pzPr * pzPi) / (momPr * momPi), -1.f, 1.f); + histos.fill(HIST("QAChecks/Pair/hOpeningAngleBefore"), std::acos(cosAngle)); + } + + void fillPairQAAfter(float pxPr, float pyPr, float pzPr, float pxPi, float pyPi, float pzPi, float momPr, float momPi) + { + const float cosAngle = std::clamp((pxPr * pxPi + pyPr * pyPi + pzPr * pzPi) / (momPr * momPi), -1.f, 1.f); + histos.fill(HIST("QAChecks/Pair/hOpeningAngleAfter"), std::acos(cosAngle)); + } + + struct TrackCandidate { + float px, py, pz, pt, eta, phi, dcaXY, dcaZ; + float mom; + int sign; + bool hasTOF; + float tpcNSigmaPr, tpcNSigmaPi; + uint64_t globalIndex; + }; + + template + void buildCandidatePools(TrackCollection const& tracks, std::vector& protons, std::vector& pions) + { + protons.clear(); + pions.clear(); + for (auto const& track : tracks) { + if (!passesBasicTrackSelection(track)) + continue; + const float mom = RecoDecay::p(track.px(), track.py(), track.pz()); + const bool isProton = passesProtonPID(track, mom) && passesProtonDCASelection(track) && + !(requireTOFForProton && !track.hasTOF()); + const bool isPion = passesPionPID(track, mom) && passesPionDCASelection(track) && + !(requireTOFForPion && !track.hasTOF()); + if (!isProton && !isPion) + continue; + TrackCandidate cand; + cand.px = track.px(); + cand.py = track.py(); + cand.pz = track.pz(); + cand.pt = track.pt(); + cand.eta = track.eta(); + cand.phi = track.phi(); + cand.dcaXY = track.dcaXY(); + cand.dcaZ = track.dcaZ(); + cand.mom = mom; + cand.sign = track.sign(); + cand.hasTOF = track.hasTOF(); + cand.tpcNSigmaPr = track.tpcNSigmaPr(); + cand.tpcNSigmaPi = track.tpcNSigmaPi(); + cand.globalIndex = track.globalIndex(); + if (isProton) + protons.push_back(cand); + if (isPion) + pions.push_back(cand); + } + } + + template + void fillInvariantMassHistogramsFromPools( + std::vector const& protonPool, + std::vector const& pionPool, + float centrality, + bool fillPairQA = false) + { + for (auto const& protonCand : protonPool) { + for (auto const& pionCand : pionPool) { + + if constexpr (!isMixed) { + if (protonCand.globalIndex == pionCand.globalIndex) + continue; + } + + const float pxPr = protonCand.px, pyPr = protonCand.py, pzPr = protonCand.pz; + const float pxPi = pionCand.px, pyPi = pionCand.py, pzPi = pionCand.pz; + + if constexpr (!isMixed) { + if (fillPairQA) + fillPairQABefore(pxPr, pyPr, pzPr, pxPi, pyPi, pzPi, protonCand.mom, pionCand.mom); + } + + if (applyDeepAngleCut) { + const float cosAngle = std::clamp( + (pxPr * pxPi + pyPr * pyPi + pzPr * pzPi) / (protonCand.mom * pionCand.mom), -1.f, 1.f); + if (std::acos(cosAngle) < static_cast(deepAngleCutValue)) + continue; + } + + const std::array, 2> bothMomenta = { + std::array{pxPr, pyPr, pzPr}, + std::array{pxPi, pyPi, pzPi}}; + const float pairMass = RecoDecay::m(bothMomenta, std::array{massProton, massPion}); + const float pairPt = RecoDecay::pt(std::array{pxPr + pxPi, pyPr + pyPi}); + const float pairY = RecoDecay::y(std::array{pxPr + pxPi, pyPr + pyPi, pzPr + pzPi}, pairMass); + + if (pairY < cfgMinY || pairY > cfgMaxY) + continue; + + if constexpr (!isMixed) { + if (fillPairQA) + fillPairQAAfter(pxPr, pyPr, pzPr, pxPi, pyPi, pzPi, protonCand.mom, pionCand.mom); + } + + if constexpr (isMixed) { + fillDeltaHistogramMixedEvent(protonCand.sign, pionCand.sign, pairPt, pairMass, centrality, pairY); + } else { + fillDeltaHistogramSameEvent(protonCand.sign, pionCand.sign, pairPt, pairMass, centrality, pairY); + if (enableRotationalBackground) + fillRotationalBackground(protonCand.sign, pionCand.sign, + pxPr, pyPr, pzPr, pionCand.phi, pionCand.pt, pzPi, centrality); + } + } + } + } + + Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; + Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPt); + + using EventCandidates = soa::Filtered>; + using TrackCandidates = soa::Filtered>; + using TrackCandidatesMC = soa::Filtered>; + + Preslice perCol = aod::track::collisionId; + Preslice perColMC = aod::track::collisionId; + + // ── Event-mixing binning policies (one per centrality estimator) ──────────── + using BinningTypeFT0M = ColumnBinningPolicy; + using BinningTypeFT0A = ColumnBinningPolicy; + using BinningTypeFT0C = ColumnBinningPolicy; + using BinningTypeFV0A = ColumnBinningPolicy; + using BinningTypeNTPV = ColumnBinningPolicy; + + BinningTypeFT0M binningFT0M{{cfgVtxAxis, cfgCentAxis}, true}; + BinningTypeFT0A binningFT0A{{cfgVtxAxis, cfgCentAxis}, true}; + BinningTypeFT0C binningFT0C{{cfgVtxAxis, cfgCentAxis}, true}; + BinningTypeFV0A binningFV0A{{cfgVtxAxis, cfgCentAxis}, true}; + BinningTypeNTPV binningNTPV{{cfgVtxAxis, cfgCentAxis}, true}; + + SameKindPair pairFT0M{binningFT0M, cfgNoMixedEvents, -1, &cache}; + SameKindPair pairFT0A{binningFT0A, cfgNoMixedEvents, -1, &cache}; + SameKindPair pairFT0C{binningFT0C, cfgNoMixedEvents, -1, &cache}; + SameKindPair pairFV0A{binningFV0A, cfgNoMixedEvents, -1, &cache}; + SameKindPair pairNTPV{binningNTPV, cfgNoMixedEvents, -1, &cache}; + + std::vector mProtonPool{}; + std::vector mPionPool{}; + + void processSameEvent(EventCandidates const& collisions, + TrackCandidates const& tracks, + aod::BCs const&) + { + for (auto const& collision : collisions) { + if (!passesEventSelection(collision)) + continue; + const float centrality = getCentrality(collision); + const int occupancy = collision.trackOccupancyInTimeRange(); + + histos.fill(HIST("OccupancyQA/hOccupancyVsCentralityBefore"), centrality, occupancy); + histos.fill(HIST("OccupancyQA/hOccupancyVsVtxZBefore"), collision.posZ(), occupancy); + histos.fill(HIST("Event/hNcontributor"), collision.numContrib()); + histos.fill(HIST("Event/hVtxZ"), collision.posZ()); + histos.fill(HIST("Event/hCentrality"), centrality); + histos.fill(HIST("Event/hOccupancy"), occupancy); + histos.fill(HIST("CentQA/hCentralityVsVtxZ"), collision.posZ(), centrality); + histos.fill(HIST("CentQA/hCentralityVsOccupancy"), occupancy, centrality); + histos.fill(HIST("CentQA/hEventCountVsCentrality"), centrality); + histos.fill(HIST("OccupancyQA/hOccupancyVsCentralityAfter"), centrality, occupancy); + histos.fill(HIST("OccupancyQA/hOccupancyVsVtxZAfter"), collision.posZ(), occupancy); + + const uint64_t collIdx = collision.globalIndex(); + auto perColTracks = tracks.sliceBy(perCol, collIdx); + perColTracks.bindExternalIndices(&tracks); + + for (auto const& track : perColTracks) { + const float mom = RecoDecay::p(track.px(), track.py(), track.pz()); + histos.fill(HIST("QAbefore/Proton/tpcNSigmaVsMomentum"), mom, track.tpcNSigmaPr()); + if (track.hasTOF()) { + histos.fill(HIST("QAbefore/Proton/tofNSigmaVsMomentum"), mom, track.tofNSigmaPr()); + histos.fill(HIST("QAbefore/Proton/tofNSigmaVsTPCNSigma"), track.tpcNSigmaPr(), track.tofNSigmaPr()); + } + histos.fill(HIST("QAbefore/Pion/tpcNSigmaVsMomentum"), mom, track.tpcNSigmaPi()); + if (track.hasTOF()) { + histos.fill(HIST("QAbefore/Pion/tofNSigmaVsMomentum"), mom, track.tofNSigmaPi()); + histos.fill(HIST("QAbefore/Pion/tofNSigmaVsTPCNSigma"), track.tpcNSigmaPi(), track.tofNSigmaPi()); + } + if (!passesBasicTrackSelection(track)) + continue; + if (passesProtonPID(track, mom) && passesProtonDCASelection(track)) + fillQAProton(track, mom, centrality); + if (passesPionPID(track, mom) && passesPionDCASelection(track)) + fillQAPion(track, mom, centrality); + } + buildCandidatePools(perColTracks, mProtonPool, mPionPool); + fillInvariantMassHistogramsFromPools(mProtonPool, mPionPool, centrality, true); + } + } + PROCESS_SWITCH(DeltaAnalysis, processSameEvent, "Process same event", true); + + // ── Shared mixed-event logic ───────────────────────────────────────────── + template + void runMixedEvent(PairType& mixingPair) + { + for (auto const& [c1, tracks1, c2, tracks2] : mixingPair) { + if (!passesEventSelection(c1) || !passesEventSelection(c2)) + continue; + const float centrality = getCentrality(c1); + std::vector protonPool1, pionPool1, protonPool2, pionPool2; + buildCandidatePools(tracks1, protonPool1, pionPool1); + buildCandidatePools(tracks2, protonPool2, pionPool2); + fillInvariantMassHistogramsFromPools(protonPool1, pionPool2, centrality); + fillInvariantMassHistogramsFromPools(protonPool2, pionPool1, centrality); + } + } + + void processMixedEvent(EventCandidates const&, TrackCandidates const&) + { + switch (cfgCentralityEstimator) { + case delta_analysis::kFT0M: + runMixedEvent(pairFT0M); + break; + case delta_analysis::kFT0A: + runMixedEvent(pairFT0A); + break; + case delta_analysis::kFT0C: + runMixedEvent(pairFT0C); + break; + case delta_analysis::kFV0A: + runMixedEvent(pairFV0A); + break; + case delta_analysis::kNTPV: + runMixedEvent(pairNTPV); + break; + default: + runMixedEvent(pairFT0M); + break; + } + } + PROCESS_SWITCH(DeltaAnalysis, processMixedEvent, "Process mixed event", true); + + void processMC(soa::Join const& collisions, aod::BCs const&, TrackCandidatesMC const& tracks, aod::McParticles const& mcParticles) + { + constexpr float kGenCentrality = 1.f; + + for (auto const& collision : collisions) { + if (!collision.sel8() || std::abs(collision.posZ()) > cfgCutVertex) + continue; + const float centrality = getCentrality(collision); + histos.fill(HIST("Event/hNcontributor"), collision.numContrib()); + histos.fill(HIST("Event/hVtxZ"), collision.posZ()); + + const uint64_t collIdx = collision.globalIndex(); + auto perColTracks = tracks.sliceBy(perColMC, collIdx); + perColTracks.bindExternalIndices(&tracks); + + for (auto const& t0 : perColTracks) { + if (!passesBasicTrackSelection(t0) || !t0.has_mcParticle()) + continue; + const auto mcTrack = t0.mcParticle(); + const float mom = RecoDecay::p(t0.px(), t0.py(), t0.pz()); + if (std::abs(mcTrack.pdgCode()) == delta_analysis::PdgProton && passesProtonPID(t0, mom)) { + for (const auto& mother : mcTrack.mothers_as()) { + if (std::abs(mother.pdgCode()) == delta_analysis::PdgDeltaPlusPlus) { + if (t0.sign() < 0) + histos.fill(HIST("QAChecks/hRecProtonAntiDeltaPlusPlus"), t0.pt()); + else + histos.fill(HIST("QAChecks/hRecProtonDeltaPlusPlus"), t0.pt()); + } else if (std::abs(mother.pdgCode()) == delta_analysis::PdgDeltaZero) { + if (t0.sign() < 0) + histos.fill(HIST("QAChecks/hRecProtonAntiDeltaZero"), t0.pt()); + else + histos.fill(HIST("QAChecks/hRecProtonDeltaZero"), t0.pt()); + } + } + } + if (std::abs(mcTrack.pdgCode()) == delta_analysis::PdgPion && passesPionPID(t0, mom)) { + for (const auto& mother : mcTrack.mothers_as()) { + if (std::abs(mother.pdgCode()) == delta_analysis::PdgDeltaPlusPlus) { + if (t0.sign() < 0) + histos.fill(HIST("QAChecks/hRecPionAntiDeltaPlusPlus"), t0.pt()); + else + histos.fill(HIST("QAChecks/hRecPionDeltaPlusPlus"), t0.pt()); + } else if (std::abs(mother.pdgCode()) == delta_analysis::PdgDeltaZero) { + if (t0.sign() < 0) + histos.fill(HIST("QAChecks/hRecPionAntiDeltaZero"), t0.pt()); + else + histos.fill(HIST("QAChecks/hRecPionDeltaZero"), t0.pt()); + } + } + } + } + + for (auto const& [t0, t1] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(perColTracks, perColTracks))) { + if (t0.globalIndex() == t1.globalIndex()) + continue; + if (!passesBasicTrackSelection(t0) || !passesBasicTrackSelection(t1)) + continue; + const float momT0 = RecoDecay::p(t0.px(), t0.py(), t0.pz()); + const float momT1 = RecoDecay::p(t1.px(), t1.py(), t1.pz()); + if (!passesProtonPID(t0, momT0) || !passesPionPID(t1, momT1)) + continue; + if (!passesProtonDCASelection(t0) || !passesPionDCASelection(t1)) + continue; + if (!t0.has_mcParticle() || !t1.has_mcParticle()) + continue; + const auto mcProton = t0.mcParticle(); + const auto mcPion = t1.mcParticle(); + if (std::abs(mcProton.pdgCode()) != delta_analysis::PdgProton) + continue; + if (std::abs(mcPion.pdgCode()) != delta_analysis::PdgPion) + continue; + bool foundMother = false; + for (const auto& motherPr : mcProton.mothers_as()) { + for (const auto& motherPi : mcPion.mothers_as()) { + if (motherPr != motherPi) + continue; + if (std::abs(motherPr.pdgCode()) != delta_analysis::PdgDeltaPlusPlus && + std::abs(motherPr.pdgCode()) != delta_analysis::PdgDeltaZero) + continue; + foundMother = true; + break; + } + if (foundMother) + break; + } + if (!foundMother) + continue; + const std::array, 2> momenta = { + std::array{t0.px(), t0.py(), t0.pz()}, + std::array{t1.px(), t1.py(), t1.pz()}}; + const float pairMass = RecoDecay::m(momenta, std::array{massProton, massPion}); + const float pairPt = RecoDecay::pt(std::array{t0.px() + t1.px(), t0.py() + t1.py()}); + const float pairY = RecoDecay::y(std::array{t0.px() + t1.px(), t0.py() + t1.py(), t0.pz() + t1.pz()}, pairMass); + if (pairY < cfgMinY || pairY > cfgMaxY) + continue; + fillDeltaHistogramMC(t0.sign(), t1.sign(), pairPt, pairMass, centrality, pairY); + } + } + + for (auto const& mcParticle : mcParticles) { + const int pdg = mcParticle.pdgCode(); + if (std::abs(pdg) != delta_analysis::PdgDeltaPlusPlus && + std::abs(pdg) != delta_analysis::PdgDeltaZero) + continue; + if (mcParticle.y() < cfgMinY || mcParticle.y() > cfgMaxY) + continue; + const auto daughters = mcParticle.daughters_as(); + bool hasPr = false, hasPi = false; + float ptPr = -999.f, ptPi = -999.f; + double eTotal = 0.; + for (const auto& d : daughters) { + if (std::abs(d.pdgCode()) == delta_analysis::PdgProton) { + hasPr = true; + ptPr = d.pt(); + eTotal += d.e(); + } else if (std::abs(d.pdgCode()) == delta_analysis::PdgPion) { + hasPi = true; + ptPi = d.pt(); + eTotal += d.e(); + } + } + if (!hasPr || !hasPi) + continue; + const float pSq = mcParticle.p() * mcParticle.p(); + const float genMass = std::sqrt(std::max(0., eTotal * eTotal - static_cast(pSq))); + const float genPt = mcParticle.pt(); + const float genY = mcParticle.y(); + if (pdg == delta_analysis::PdgDeltaPlusPlus) { + histos.fill(HIST("Analysis/hDeltaPlusPlusInvMassGen"), genPt, genMass); + histos.fill(HIST("THnSparse/hDeltaPlusPlusGen"), genPt, kGenCentrality, genY); + histos.fill(HIST("QAChecks/hGenProtonDeltaPlusPlus"), ptPr); + histos.fill(HIST("QAChecks/hGenPionDeltaPlusPlus"), ptPi); + } else if (pdg == -delta_analysis::PdgDeltaPlusPlus) { + histos.fill(HIST("Analysis/hAntiDeltaPlusPlusInvMassGen"), genPt, genMass); + histos.fill(HIST("THnSparse/hAntiDeltaPlusPlusGen"), genPt, kGenCentrality, genY); + histos.fill(HIST("QAChecks/hGenProtonAntiDeltaPlusPlus"), ptPr); + histos.fill(HIST("QAChecks/hGenPionAntiDeltaPlusPlus"), ptPi); + } else if (pdg == delta_analysis::PdgDeltaZero) { + histos.fill(HIST("Analysis/hDeltaZeroInvMassGen"), genPt, genMass); + histos.fill(HIST("THnSparse/hDeltaZeroGen"), genPt, kGenCentrality, genY); + histos.fill(HIST("QAChecks/hGenProtonDeltaZero"), ptPr); + histos.fill(HIST("QAChecks/hGenPionDeltaZero"), ptPi); + } else if (pdg == -delta_analysis::PdgDeltaZero) { + histos.fill(HIST("Analysis/hAntiDeltaZeroInvMassGen"), genPt, genMass); + histos.fill(HIST("THnSparse/hAntiDeltaZeroGen"), genPt, kGenCentrality, genY); + histos.fill(HIST("QAChecks/hGenProtonAntiDeltaZero"), ptPr); + histos.fill(HIST("QAChecks/hGenPionAntiDeltaZero"), ptPi); + } + } + } + PROCESS_SWITCH(DeltaAnalysis, processMC, "Process MC", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + // FIX name/o2-task: TaskName is redundant when it equals the derived device name + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGLF/Tasks/Resonances/deltaanalysis.cxx b/PWGLF/Tasks/Resonances/deltaanalysis.cxx deleted file mode 100644 index a121980f67e..00000000000 --- a/PWGLF/Tasks/Resonances/deltaanalysis.cxx +++ /dev/null @@ -1,613 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. -// Analysis task for delta analysis - -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/PIDResponseTOF.h" -#include "Common/DataModel/PIDResponseTPC.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include -#include -#include - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; - -namespace -{ -constexpr float pionMass = o2::constants::physics::MassPionCharged; -constexpr float protonMass = o2::constants::physics::MassProton; -constexpr int deltaPlusPlusPDG = 2224; -constexpr int deltaZeroPDG = 2114; -constexpr int protonPDG = 2212; -constexpr int pionPDG = 211; -} // namespace - -struct deltaAnalysis { - - SliceCache cache; - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - // events - Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; - // track - Configurable cfgCutPt{"cfgCutPt", 0.2, "Pt cut on daughter track"}; - Configurable cfgCutEta{"cfgCutEta", 0.8, "Eta cut on daughter track"}; - Configurable cfgCutY{"cfgCutY", 0.5, "Y cut on reconstructed delta"}; - Configurable cfgCutDCAxy{"cfgCutDCAxy", 2.0f, "DCAxy range for tracks"}; - Configurable cfgCutDCAz{"cfgCutDCAz", 2.0f, "DCAz range for tracks"}; - Configurable nsigmaCutTPC{"nsigmacutTPC", 3.0, "Value of the TPC Nsigma cut"}; - Configurable nsigmaCutTOF{"nsigmaCutTOF", 3.0, "Value of the TOF Nsigma cut"}; - Configurable cfgNoMixedEvents{"cfgNoMixedEvents", 5, "Number of mixed events per event"}; - Configurable cfgCutPtProtonTPC{"cfgCutPtProtonTPC", 0.8, "Pt cut of proton in TPC"}; - Configurable cfgCutPtPionTPC{"cfgCutPtPionTPC", 0.8, "Pt cut of pion in TPC"}; - - // Histogram axes - AxisSpec nSigmaTPCaxis = {100, -5., 5., "n#sigma_{TPC}"}; - AxisSpec nSigmaTOFaxis = {100, -5., 5., "n#sigma_{TOF}"}; - ConfigurableAxis cfgPtAxis{"cfgPtAxis", {VARIABLE_WIDTH, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4.}, "#it{p}_{T} (GeV/#it{c})"}; - ConfigurableAxis cfgMassAxis{"cfgDeltaPlusPlusAxis", {75, 1.05, 1.8}, ""}; - - void init(o2::framework::InitContext&) - { - // Define axes - AxisSpec deltaPlusPlusAxis{cfgMassAxis, "m(p + #pi^{+}) (GeV/#it{c}^{2})"}; - AxisSpec antiDeltaPlusPlusAxis{cfgMassAxis, "m(#bar{p} + #pi^{-}) (GeV/#it{c}^{2})"}; - AxisSpec deltaZeroAxis{cfgMassAxis, "m(p + #pi^{-}) (GeV/#it{c}^{2})"}; - AxisSpec antiDeltaZeroAxis{cfgMassAxis, "m(#bar{p} + #pi^{+}) (GeV/#it{c}^{2})"}; - AxisSpec ptAxis{cfgPtAxis, "#it{p}_{T} (GeV/#it{c})"}; - - // Collision - histos.add("hCentrality", "Centrality distribution", kTH1F, {{2001, -0.5, 2000.5}}); - histos.add("hVtxZ", "Vertex distribution in Z;Z (cm)", kTH1F, {{400, -20.0, 20.0}}); - histos.add("hNcontributor", "Number of primary vertex contributors; n", kTH1F, {{2001, -0.5f, 2000.5f}}); - - // Single track - histos.add("hPiPlusDCAxy", "DCA_{xy} distribution for #pi^{+}; DCA_{xy} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hPiPlusDCAz", "DCA_{z} distribution for #pi^{+}; DCA_{z} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hPiPlusNsigmaTPCvsPt", "n#sigma_{TPC} distribution vs #it{p}_{T} for #pi^{+}", kTH2F, {ptAxis, nSigmaTPCaxis}); - histos.add("hPiPlusNsigmaTPCvsPt_TPC_only", "n#sigma_{TPC} distribution vs #it{p}_{T} for #pi^{+}", kTH2F, {ptAxis, nSigmaTPCaxis}); - histos.add("hPiPlusNsigmaTOFvsPt", "n#sigma_{TOF} distribution vs #it{p}_{T} for #pi^{+}", kTH2F, {ptAxis, nSigmaTOFaxis}); - - histos.add("hPiMinusDCAxy", "DCA_{xy} distribution for #pi^{-}; DCA_{xy} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hPiMinusDCAz", "DCA_{z} distribution for #pi^{-}; DCA_{z} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hPiMinusNsigmaTPCvsPt", "n#sigma_{TPC} distribution vs #it{p}_{T} for #pi^{-}", kTH2F, {ptAxis, nSigmaTPCaxis}); - histos.add("hPiMinusNsigmaTPCvsPt_TPC_only", "n#sigma_{TPC} distribution vs #it{p}_{T} for #pi^{-}", kTH2F, {ptAxis, nSigmaTPCaxis}); - histos.add("hPiMinusNsigmaTOFvsPt", "n#sigma_{TOF} distribution vs #it{p}_{T} for #pi^{-}", kTH2F, {ptAxis, nSigmaTOFaxis}); - - histos.add("hPrPlusDCAxy", "DCA_{xy} distribution for p; DCA_{xy} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hPrPlusDCAz", "DCA_{z} distribution for p; DCA_{z} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hPrPlusNsigmaTPCvsPt", "n#sigma_{TPC} distribution vs #it{p}_{T} for p", kTH2F, {ptAxis, nSigmaTPCaxis}); - histos.add("hPrPlusNsigmaTPCvsPt_TPC_only", "n#sigma_{TPC} distribution vs #it{p}_{T} for p", kTH2F, {ptAxis, nSigmaTPCaxis}); - histos.add("hPrPlusNsigmaTOFvsPt", "n#sigma_{TOF} distribution vs #it{p}_{T} for p", kTH2F, {ptAxis, nSigmaTOFaxis}); - - histos.add("hPrMinusDCAxy", "DCA_{xy} distribution for #bar{p}; DCA_{xy} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hPrMinusDCAz", "DCA_{z} distribution for #bar{p}; DCA_{z} (cm)", kTH1F, {{200, -1.0f, 1.0f}}); - histos.add("hPrMinusNsigmaTPCvsPt", "n#sigma_{TPC} distribution vs #it{p}_{T} for #bar{p}", kTH2F, {ptAxis, nSigmaTPCaxis}); - histos.add("hPrMinusNsigmaTPCvsPt_TPC_only", "n#sigma_{TPC} distribution vs #it{p}_{T} for #bar{p}", kTH2F, {ptAxis, nSigmaTPCaxis}); - histos.add("hPrMinusNsigmaTOFvsPt", "n#sigma_{TOF} distribution vs #it{p}_{T} for #bar{p}", kTH2F, {ptAxis, nSigmaTOFaxis}); - - // Deltas - histos.add("hDeltaPlusPlusInvMass", "Invariant mass distribution for #Delta^{++}", kTH2F, {ptAxis, deltaPlusPlusAxis}); - histos.add("hAntiDeltaPlusPlusInvMass", "Invariant mass distribution for #bar{#Delta^{++}}", kTH2F, {ptAxis, antiDeltaPlusPlusAxis}); - - histos.add("hDeltaZeroInvMass", "Invariant mass distribution for #Delta^{0}", kTH2F, {ptAxis, deltaZeroAxis}); - histos.add("hAntiDeltaZeroInvMass", "Invariant mass distribution for #bar{#Delta^{0}}", kTH2F, {ptAxis, antiDeltaZeroAxis}); - - if (doprocessMixedEvent) { - // Deltas - Event mixing - histos.add("hDeltaPlusPlusInvMassEM", "Invariant mass distribution for #Delta^{++} - event mixing", kTH2F, {ptAxis, deltaPlusPlusAxis}); - histos.add("hAntiDeltaPlusPlusInvMassEM", "Invariant mass distribution for #bar{#Delta^{++}} - event mixing", kTH2F, {ptAxis, antiDeltaPlusPlusAxis}); - - histos.add("hDeltaZeroInvMassEM", "Invariant mass distribution for #Delta^{0} - event mixing", kTH2F, {ptAxis, deltaZeroAxis}); - histos.add("hAntiDeltaZeroInvMassEM", "Invariant mass distribution for #bar{#Delta^{0}} - event mixing", kTH2F, {ptAxis, antiDeltaZeroAxis}); - } - - if (doprocessMC) { - // generated quantities - histos.add("hDeltaPlusPlusInvMassGen", "Invariant mass distribution for #Delta^{++} - generated", kTH2F, {ptAxis, deltaPlusPlusAxis}); - histos.add("hAntiDeltaPlusPlusInvMassGen", "Invariant mass distribution for #bar{#Delta^{++}} - generated", kTH2F, {ptAxis, antiDeltaPlusPlusAxis}); - - histos.add("hDeltaZeroInvMassGen", "Invariant mass distribution for #Delta^{0} - generated", kTH2F, {ptAxis, deltaZeroAxis}); - histos.add("hAntiDeltaZeroInvMassGen", "Invariant mass distribution for #bar{#Delta^{0}} - generated", kTH2F, {ptAxis, antiDeltaZeroAxis}); - - histos.add("hRecProtonAntiDeltaPlusPlus", "Proton from #bar{#Delta^{++}}", kTH1F, {ptAxis}); - histos.add("hRecProtonDeltaPlusPlus", "Proton from #Delta^{++}", kTH1F, {ptAxis}); - histos.add("hRecProtonAntiDeltaZero", "Proton from #bar{#Delta^{0}}", kTH1F, {ptAxis}); - histos.add("hRecProtonDeltaZero", "Proton from #Delta^{0}", kTH1F, {ptAxis}); - histos.add("hRecPionAntiDeltaPlusPlus", "Pion from #bar{#Delta^{++}}", kTH1F, {ptAxis}); - histos.add("hRecPionDeltaPlusPlus", "Pion from #Delta^{++}", kTH1F, {ptAxis}); - histos.add("hRecPionAntiDeltaZero", "Pion from #bar{#Delta^{0}}", kTH1F, {ptAxis}); - histos.add("hRecPionDeltaZero", "Pion from #Delta^{0}", kTH1F, {ptAxis}); - - histos.add("hGenProtonAntiDeltaPlusPlus", "Proton from #bar{#Delta^{++}}", kTH1F, {ptAxis}); - histos.add("hGenProtonDeltaPlusPlus", "Proton from #Delta^{++}", kTH1F, {ptAxis}); - histos.add("hGenProtonAntiDeltaZero", "Proton from #bar{#Delta^{0}}", kTH1F, {ptAxis}); - histos.add("hGenProtonDeltaZero", "Proton from #Delta^{0}", kTH1F, {ptAxis}); - histos.add("hGenPionAntiDeltaPlusPlus", "Pion from #bar{#Delta^{++}}", kTH1F, {ptAxis}); - histos.add("hGenPionDeltaPlusPlus", "Pion from #Delta^{++}", kTH1F, {ptAxis}); - histos.add("hGenPionAntiDeltaZero", "Pion from #bar{#Delta^{0}}", kTH1F, {ptAxis}); - histos.add("hGenPionDeltaZero", "Pion from #Delta^{0}", kTH1F, {ptAxis}); - } - } - - template - bool selectionTrack(const T& track) - { - // if (!track.isGlobalTrack()) { - // return false; - // } - if (track.itsNCls() < 5) { - return false; - } - if (track.tpcNClsShared() > 0) { - return false; - } - if (track.tpcNClsFound() < 70) { - return false; - } - if (track.tpcNClsCrossedRows() < 70) { - return false; - } - if (track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable()) { - return false; - } - if (track.tpcChi2NCl() > 4.f) { - return false; - } - if (track.itsChi2NCl() > 36.f) { - return false; - } - return true; - } - - template - bool selectionPIDPion(const T& track) - { - if (track.hasTOF()) { - if (std::abs(track.tofNSigmaPi()) < nsigmaCutTOF && std::abs(track.tpcNSigmaPi()) < nsigmaCutTPC) { - if (track.sign() > 0) { - histos.fill(HIST("hPiPlusNsigmaTPCvsPt"), track.pt(), track.tpcNSigmaPi()); - histos.fill(HIST("hPiPlusNsigmaTOFvsPt"), track.pt(), track.tofNSigmaPi()); - histos.fill(HIST("hPiPlusDCAxy"), track.dcaXY()); - histos.fill(HIST("hPiPlusDCAz"), track.dcaZ()); - } else { - histos.fill(HIST("hPiMinusNsigmaTPCvsPt"), track.pt(), track.tpcNSigmaPi()); - histos.fill(HIST("hPiMinusNsigmaTOFvsPt"), track.pt(), track.tofNSigmaPi()); - histos.fill(HIST("hPiMinusDCAxy"), track.dcaXY()); - histos.fill(HIST("hPiMinusDCAz"), track.dcaZ()); - } - return true; - } - } else if (std::abs(track.tpcNSigmaPi()) < nsigmaCutTPC) { - if (track.sign() > 0) { - histos.fill(HIST("hPiPlusNsigmaTPCvsPt"), track.pt(), track.tpcNSigmaPi()); - histos.fill(HIST("hPiPlusNsigmaTPCvsPt_TPC_only"), track.pt(), track.tpcNSigmaPi()); - histos.fill(HIST("hPiPlusDCAxy"), track.dcaXY()); - histos.fill(HIST("hPiPlusDCAz"), track.dcaZ()); - } else { - histos.fill(HIST("hPiMinusNsigmaTPCvsPt"), track.pt(), track.tpcNSigmaPi()); - histos.fill(HIST("hPiMinusNsigmaTPCvsPt_TPC_only"), track.pt(), track.tpcNSigmaPi()); - histos.fill(HIST("hPiMinusDCAxy"), track.dcaXY()); - histos.fill(HIST("hPiMinusDCAz"), track.dcaZ()); - } - if (track.pt() < cfgCutPtPionTPC) { - return true; - } - } - return false; - } - - template - bool selectionPIDProton(const T& track) - { - if (track.hasTOF()) { - if (std::abs(track.tofNSigmaPr()) < nsigmaCutTOF && std::abs(track.tpcNSigmaPr()) < nsigmaCutTPC) { - if (track.sign() > 0) { - histos.fill(HIST("hPrPlusNsigmaTPCvsPt"), track.pt(), track.tpcNSigmaPr()); - histos.fill(HIST("hPrPlusNsigmaTOFvsPt"), track.pt(), track.tofNSigmaPr()); - histos.fill(HIST("hPrPlusDCAxy"), track.dcaXY()); - histos.fill(HIST("hPrPlusDCAz"), track.dcaZ()); - } else { - histos.fill(HIST("hPrMinusNsigmaTPCvsPt"), track.pt(), track.tpcNSigmaPr()); - histos.fill(HIST("hPrMinusNsigmaTOFvsPt"), track.pt(), track.tofNSigmaPr()); - histos.fill(HIST("hPrMinusDCAxy"), track.dcaXY()); - histos.fill(HIST("hPrMinusDCAz"), track.dcaZ()); - } - return true; - } - } else if (std::abs(track.tpcNSigmaPr()) < nsigmaCutTPC) { - if (track.sign() > 0) { - histos.fill(HIST("hPrPlusNsigmaTPCvsPt"), track.pt(), track.tpcNSigmaPr()); - histos.fill(HIST("hPrPlusNsigmaTPCvsPt_TPC_only"), track.pt(), track.tpcNSigmaPr()); - histos.fill(HIST("hPrPlusDCAxy"), track.dcaXY()); - histos.fill(HIST("hPrPlusDCAz"), track.dcaZ()); - } else { - histos.fill(HIST("hPrMinusNsigmaTPCvsPt"), track.pt(), track.tpcNSigmaPr()); - histos.fill(HIST("hPrMinusNsigmaTPCvsPt_TPC_only"), track.pt(), track.tpcNSigmaPr()); - histos.fill(HIST("hPrMinusDCAxy"), track.dcaXY()); - histos.fill(HIST("hPrMinusDCAz"), track.dcaZ()); - } - if (track.pt() < cfgCutPtProtonTPC) { - return true; - } - } - return false; - } - - Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; - Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPt); - Filter DCAcutFilter = (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz); - - using EventCandidates = soa::Filtered>; - - using TrackCandidates = soa::Filtered>; - - using TrackCandidatesMC = soa::Filtered>; - - Preslice perCol = aod::track::collisionId; - Preslice perColMC = aod::track::collisionId; - - // Binning for event mixing - ConfigurableAxis cfgCentAxis{"cfgCentAxis", {VARIABLE_WIDTH, 0., 0.01, 0.1, 1.0, 5.0, 10., 15., 20., 30., 40., 50., 70., 100.0, 105.}, "Binning of the centrality axis"}; - ConfigurableAxis cfgVtxAxis{"cfgVtxAxis", {VARIABLE_WIDTH, -20, -15, -10, -7, -5, -3, -2, -1, 0, 1, 2, 3, 5, 7, 10, 15, 20}, "Mixing bins - z-vertex"}; - - using BinningType = ColumnBinningPolicy; - - BinningType binningOnPositions{{cfgVtxAxis}, true}; - SameKindPair pair{binningOnPositions, cfgNoMixedEvents, -1, &cache}; - - void processSameEvent(EventCandidates const& collisions, TrackCandidates const& tracks, aod::BCs const&) - { - for (auto& collision : collisions) { - if (!collision.sel8()) { - continue; - } - histos.fill(HIST("hNcontributor"), collision.numContrib()); - histos.fill(HIST("hVtxZ"), collision.posZ()); - - const uint64_t collIdx = collision.globalIndex(); - auto perColTracks = tracks.sliceBy(perCol, collIdx); - perColTracks.bindExternalIndices(&tracks); - - for (auto& [t0, t1] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(perColTracks, perColTracks))) { - - if (t0.globalIndex() == t1.globalIndex()) { - continue; - } - if (!selectionTrack(t0) || !selectionTrack(t1)) { - continue; - } - if (selectionPIDProton(t0)) { - - TLorentzVector proton; - proton.SetPtEtaPhiM(t0.pt(), t0.eta(), t0.phi(), protonMass); - - if (selectionPIDPion(t1)) { - - TLorentzVector pion; - pion.SetPtEtaPhiM(t1.pt(), t1.eta(), t1.phi(), pionMass); - - TLorentzVector delta = proton + pion; - - if (abs(delta.Rapidity()) > cfgCutY) { - continue; - } - - if (t0.sign() < 0) { - if (t1.sign() < 0) { - histos.fill(HIST("hAntiDeltaPlusPlusInvMass"), delta.Pt(), delta.M()); - } else { - histos.fill(HIST("hAntiDeltaZeroInvMass"), delta.Pt(), delta.M()); - } - } else { - if (t1.sign() < 0) { - histos.fill(HIST("hDeltaZeroInvMass"), delta.Pt(), delta.M()); - } else { - histos.fill(HIST("hDeltaPlusPlusInvMass"), delta.Pt(), delta.M()); - } - } - } - } - } - } - } - PROCESS_SWITCH(deltaAnalysis, processSameEvent, "Process same event", false); - - void processMixedEvent(EventCandidates const& /*collisions*/, TrackCandidates const& /*tracks*/) - { - for (auto& [c1, tracks1, c2, tracks2] : pair) { - if (!c1.sel8()) { - continue; - } - if (!c2.sel8()) { - continue; - } - - for (auto& [t0, t1] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - - if (!selectionTrack(t0) || !selectionTrack(t1)) { - continue; - } - - if (selectionPIDProton(t0)) { - - TLorentzVector proton; - proton.SetPtEtaPhiM(t0.pt(), t0.eta(), t0.phi(), protonMass); - - if (selectionPIDPion(t1)) { - - TLorentzVector pion; - pion.SetPtEtaPhiM(t1.pt(), t1.eta(), t1.phi(), pionMass); - - TLorentzVector delta = proton + pion; - - if (abs(delta.Rapidity()) > cfgCutY) { - continue; - } - - if (t0.sign() < 0) { - if (t1.sign() < 0) { - histos.fill(HIST("hAntiDeltaPlusPlusInvMassEM"), delta.Pt(), delta.M()); - } else { - histos.fill(HIST("hAntiDeltaZeroInvMassEM"), delta.Pt(), delta.M()); - } - } else { - if (t1.sign() < 0) { - histos.fill(HIST("hDeltaZeroInvMassEM"), delta.Pt(), delta.M()); - } else { - histos.fill(HIST("hDeltaPlusPlusInvMassEM"), delta.Pt(), delta.M()); - } - } - } - } - } - } - } - PROCESS_SWITCH(deltaAnalysis, processMixedEvent, "Process Mixed event", false); - - void processMC(soa::Join const& collisions, aod::BCs const&, TrackCandidatesMC const& tracks, aod::McParticles const& mcParticles) - { - - for (auto& collision : collisions) { - - if (!collision.sel8() || std::abs(collision.posZ()) > cfgCutVertex) { - continue; - } - - histos.fill(HIST("hCentrality"), 1); - histos.fill(HIST("hNcontributor"), collision.numContrib()); - histos.fill(HIST("hVtxZ"), collision.posZ()); - - const uint64_t collIdx = collision.globalIndex(); - auto perColTracks = tracks.sliceBy(perColMC, collIdx); - perColTracks.bindExternalIndices(&tracks); - - for (auto& t0 : perColTracks) { - - if (!selectionTrack(t0)) { - continue; - } - if (!t0.has_mcParticle()) { - continue; - } - - const auto mcTrack = t0.mcParticle(); - - if (std::abs(mcTrack.pdgCode()) == protonPDG) { - if (selectionPIDProton(t0)) { - for (auto& motherTrackProton : mcTrack.mothers_as()) { - if (std::abs(motherTrackProton.pdgCode()) == deltaPlusPlusPDG) { - if (t0.sign() < 0) { - histos.fill(HIST("hRecProtonAntiDeltaPlusPlus"), t0.pt()); - } else { - histos.fill(HIST("hRecProtonDeltaPlusPlus"), t0.pt()); - } - } else if (std::abs(motherTrackProton.pdgCode()) == deltaZeroPDG) { - if (t0.sign() < 0) { - histos.fill(HIST("hRecProtonAntiDeltaZero"), t0.pt()); - } else { - histos.fill(HIST("hRecProtonDeltaZero"), t0.pt()); - } - } - } - } - } else if (std::abs(mcTrack.pdgCode()) == pionPDG) { - if (selectionPIDPion(t0)) { - for (auto& motherTrackPion : mcTrack.mothers_as()) { - if (std::abs(motherTrackPion.pdgCode()) == deltaPlusPlusPDG) { - if (t0.sign() < 0) { - histos.fill(HIST("hRecPionAntiDeltaPlusPlus"), t0.pt()); - } else { - histos.fill(HIST("hRecPionDeltaPlusPlus"), t0.pt()); - } - } else if (std::abs(motherTrackPion.pdgCode()) == deltaZeroPDG) { - if (t0.sign() < 0) { - histos.fill(HIST("hRecPionDeltaZero"), t0.pt()); - } else { - histos.fill(HIST("hRecPionAntiDeltaZero"), t0.pt()); - } - } - } - } - } - } - - for (auto& [t0, t1] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(perColTracks, perColTracks))) { - - if (t0.globalIndex() == t1.globalIndex()) { - continue; - } - if (!selectionTrack(t0) || !selectionTrack(t1)) { - continue; - } - if (selectionPIDProton(t0)) { - - if (!t0.has_mcParticle()) { - continue; - } - TLorentzVector proton; - proton.SetPtEtaPhiM(t0.pt(), t0.eta(), t0.phi(), protonMass); - - if (selectionPIDPion(t1)) { - - if (!t1.has_mcParticle()) { - continue; - } - TLorentzVector pion; - pion.SetPtEtaPhiM(t1.pt(), t1.eta(), t1.phi(), pionMass); - - // select real protons and pion - const auto mcTrackProton = t0.mcParticle(); - const auto mcTrackPion = t1.mcParticle(); - - if (std::abs(mcTrackProton.pdgCode()) != protonPDG) { - continue; - } - - if (std::abs(mcTrackPion.pdgCode()) != pionPDG) { - continue; - } - - bool found_mother = false; - - for (auto& motherTrackProton : mcTrackProton.mothers_as()) { - for (auto& motherTrackPion : mcTrackPion.mothers_as()) { - if (motherTrackProton != motherTrackPion) { - continue; - } - if (std::abs(motherTrackProton.pdgCode()) != deltaPlusPlusPDG && std::abs(motherTrackProton.pdgCode()) != deltaZeroPDG) { - continue; - } - found_mother = true; - break; - } - if (found_mother) { - break; - } - } - - if (!found_mother) { - continue; - } - - TLorentzVector delta = proton + pion; - - if (abs(delta.Rapidity()) > cfgCutY) { - continue; - } - - if (t0.sign() < 0) { - if (t1.sign() < 0) { - histos.fill(HIST("hAntiDeltaPlusPlusInvMass"), delta.Pt(), delta.M()); - } else { - histos.fill(HIST("hAntiDeltaZeroInvMass"), delta.Pt(), delta.M()); - } - } else { - if (t1.sign() < 0) { - histos.fill(HIST("hDeltaZeroInvMass"), delta.Pt(), delta.M()); - } else { - histos.fill(HIST("hDeltaPlusPlusInvMass"), delta.Pt(), delta.M()); - } - } - } - } - } - } - - for (auto& mcParticle : mcParticles) { - - int pdg = mcParticle.pdgCode(); - if (std::abs(pdg) != deltaPlusPlusPDG && std::abs(pdg) != deltaZeroPDG) { - continue; - } - - if (std::abs(mcParticle.y()) > cfgCutY) { - continue; - } - - // if (!mcParticle.isPhysicalPrimary()) { - // continue; - // } - - auto daughters = mcParticle.daughters_as(); - bool daughtPr = false; - bool daughtPi = false; - float ptPr = -999.; - float ptPi = -999.; - double eDelta = 0.; - for (auto daught : daughters) { - if (std::abs(daught.pdgCode()) == protonPDG) { - daughtPr = true; - ptPr = daught.pt(); - eDelta += daught.e(); - } else if (std::abs(daught.pdgCode()) == pionPDG) { - daughtPi = true; - ptPi = daught.pt(); - eDelta += daught.e(); - } - } - - if (!daughtPr || !daughtPi) { - continue; - } - - float gen_mass = TMath::Sqrt(eDelta * eDelta - mcParticle.p() * mcParticle.p()); - - if (pdg > 0) { - if (std::abs(pdg) == deltaPlusPlusPDG) { - histos.fill(HIST("hDeltaPlusPlusInvMassGen"), mcParticle.pt(), gen_mass); - histos.fill(HIST("hGenProtonDeltaPlusPlus"), ptPr); - histos.fill(HIST("hGenPionDeltaPlusPlus"), ptPi); - } else if (std::abs(pdg) == deltaZeroPDG) { - histos.fill(HIST("hDeltaZeroInvMassGen"), mcParticle.pt(), gen_mass); - histos.fill(HIST("hGenProtonDeltaZero"), ptPr); - histos.fill(HIST("hGenPionDeltaZero"), ptPi); - } - } else { - if (std::abs(pdg) == deltaPlusPlusPDG) { - histos.fill(HIST("hAntiDeltaPlusPlusInvMassGen"), mcParticle.pt(), gen_mass); - histos.fill(HIST("hGenProtonAntiDeltaPlusPlus"), ptPr); - histos.fill(HIST("hGenPionAntiDeltaPlusPlus"), ptPi); - } else if (std::abs(pdg) == deltaZeroPDG) { - histos.fill(HIST("hAntiDeltaZeroInvMassGen"), mcParticle.pt(), gen_mass); - histos.fill(HIST("hGenProtonAntiDeltaZero"), ptPr); - histos.fill(HIST("hGenPionAntiDeltaZero"), ptPi); - } - } - } - } - PROCESS_SWITCH(deltaAnalysis, processMC, "Process MC", true); -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc, TaskName{"deltaAnalysis"})}; -}