From 16a14a07265fba73e75022a039d8a92240c942ac Mon Sep 17 00:00:00 2001 From: Veronika Barbasova Date: Mon, 8 Sep 2025 11:13:47 +0200 Subject: [PATCH] Update Signed-off-by: Veronika Barbasova --- .../Tasks/Resonances/phianalysisTHnSparse.cxx | 855 ++++++++---------- 1 file changed, 369 insertions(+), 486 deletions(-) diff --git a/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx b/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx index f9fa43c1d7f..5d16db9af39 100644 --- a/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx +++ b/PWGLF/Tasks/Resonances/phianalysisTHnSparse.cxx @@ -13,22 +13,24 @@ /// \brief Analysis of phi resonance using THnSparse histograms. /// \author Veronika Barbasova (veronika.barbasova@cern.ch) -#include -#include -#include +#include "PWGLF/Utils/rsnOutput.h" +#include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/Centrality.h" #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/TrackSelectionTables.h" + #include "Framework/ASoAHelpers.h" #include "Framework/AnalysisTask.h" #include "Framework/O2DatabasePDGPlugin.h" #include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/PID.h" -#include "PWGLF/Utils/rsnOutput.h" -// #include "TDatabasePDG.h" + +#include + +#include +#include using namespace o2; using namespace o2::analysis; @@ -46,31 +48,34 @@ struct PhianalysisTHnSparse { Configurable eventMixing{"eventMixing", "none", "Produce Event Mixing histograms of type."}; } produce; - struct : ConfigurableGroup { - Configurable verboselevel{"verboselevel", 0, "Verbose level"}; - Configurable refresh{"refresh", 0, "Freqency of print event information."}; - Configurable refreshIndex{"refreshIndex", 0, "Freqency of print event information index."}; - } verbose; - - Configurable dautherPos{"dautherPos", 3, "Particle type of the positive dauther according to ReconstructionDataFormats/PID.h (Default = Kaon)"}; - Configurable dautherNeg{"dautherNeg", 3, "Particle type of the negative dauther according to ReconstructionDataFormats/PID.h (Default = Kaon)"}; + Configurable daughterPos{"daughterPos", 3, "Particle type of the positive dauther according to ReconstructionDataFormats/PID.h (Default = Kaon)"}; + Configurable daughterNeg{"daughterNeg", 3, "Particle type of the negative dauther according to ReconstructionDataFormats/PID.h (Default = Kaon)"}; Configurable motherPDG{"motherPDG", 333, "PDG code of mother particle."}; - Configurable dautherPosPDG{"dautherPosPDG", 321, "PDG code of positive dauther particle."}; - Configurable dautherNegPDG{"dautherNegPDG", 321, "PDG code of negative dauther particle."}; + Configurable daughterPosPDG{"daughterPosPDG", 321, "PDG code of positive dauther particle."}; + Configurable daughterNegPDG{"daughterNegPDG", 321, "PDG code of negative dauther particle."}; struct : ConfigurableGroup { Configurable tpcnSigmaPos{"tpcnSigmaPos", 3.0f, "TPC NSigma cut of the positive particle."}; Configurable tpcnSigmaNeg{"tpcnSigmaNeg", 3.0f, "TPC NSigma cut of the negative particle."}; - Configurable vZ{"vZ", 10.0f, "Z vertex range."}; - Configurable y{"y", 0.5, "Rapidity cut (maximum)."}; - Configurable etatrack{"etatrack", 0.8, "Eta cut for track."}; - Configurable etapair{"etapair", 0.8, "Eta cut for pair."}; + Configurable tpcPidOnly{"tpcPidOnly", false, "Use TPC only for PID."}; + Configurable combinedNSigma{"combinedNSigma", 3.0f, "Combined NSigma cut for combined TPC and TOF NSigma cut."}; + Configurable ptTOFThreshold{"ptTOFThreshold", 0.5f, "Threshold for applying TOF."}; + Configurable rapidity{"rapidity", 0.5f, "Rapidity cut (maximum)."}; + Configurable etatrack{"etatrack", 0.8f, "Eta cut for track."}; + Configurable etapair{"etapair", 0.5f, "Eta cut for pair."}; Configurable pt{"pt", 0.15f, "Cut: Minimal value of tracks pt."}; Configurable dcaXY{"dcaXY", 1.0f, "Cut: Maximal value of tracks DCA XY."}; Configurable dcaZ{"dcaZ", 1.0f, "Cut: Maximal value of tracks DCA Z."}; + Configurable globalTrack{"globalTrack", false, "Use global track selection."}; Configurable tpcNClsFound{"tpcNClsFound", 70, "Cut: Minimal value of found TPC clasters"}; } cut; + struct : ConfigurableGroup { + Configurable verboselevel{"verboselevel", 0, "Verbose level"}; + Configurable refresh{"refresh", 0, "Freqency of print event information."}; + Configurable refreshIndex{"refreshIndex", 0, "Freqency of print event information index."}; + } verbose; + Configurable> sparseAxes{"sparseAxes", std::vector{o2::analysis::rsn::pair_axis::names}, "Axes."}; Configurable> sysAxes{"sysAxes", std::vector{o2::analysis::rsn::systematic_axis::names}, "Axes."}; @@ -80,9 +85,9 @@ struct PhianalysisTHnSparse { ConfigurableAxis multiplicityaxis{"multiplicityaxis", {50, 0., 5000.}, "Multiplicity axis binning."}; ConfigurableAxis centralityaxis{"centralityaxis", {20, 0., 100.}, "Centrality axis binning."}; ConfigurableAxis etaaxis{"etaaxis", {16., -1.0 * static_cast(cut.etatrack), static_cast(cut.etatrack)}, "Pseudorapidity axis binning."}; - ConfigurableAxis rapidityaxis{"rapidityaxis", {10., -1.0 * static_cast(cut.y), static_cast(cut.y)}, "Rapidity axis binning."}; - ConfigurableAxis nsigmaaxisPos{"nsigmaaxisPos", {1, 0., static_cast(cut.tpcnSigmaPos)}, "NSigma of positive particle axis binning in THnSparse."}; - ConfigurableAxis nsigmaaxisNeg{"nsigmaaxisNeg", {1, 0., static_cast(cut.tpcnSigmaNeg)}, "NSigma of negative particle axis binning in THnSparse."}; + ConfigurableAxis rapidityaxis{"rapidityaxis", {10., -1.0 * static_cast(cut.rapidity), static_cast(cut.rapidity)}, "Rapidity axis binning."}; + ConfigurableAxis nsigmaaxisPos{"nsigmaaxisPos", {1, -static_cast(cut.tpcnSigmaPos), static_cast(cut.tpcnSigmaPos)}, "NSigma of positive particle axis binning in THnSparse."}; + ConfigurableAxis nsigmaaxisNeg{"nsigmaaxisNeg", {1, -static_cast(cut.tpcnSigmaNeg), static_cast(cut.tpcnSigmaNeg)}, "NSigma of negative particle axis binning in THnSparse."}; // mixing using BinningTypeVzMu = ColumnBinningPolicy>; @@ -92,14 +97,10 @@ struct PhianalysisTHnSparse { ConfigurableAxis axisMultiplicityMixing{"axisMultiplicityMixing", {5, 0, 5000}, "TPC multiplicity for bin"}; ConfigurableAxis axisCentralityMixing{"axisCentralityMixing", {10, 0, 100}, "Multiplicity percentil binning for mixing"}; - // defined in DataFormats/Reconstruction/include/ReconstructionDataFormats/PID.h - float massPos = o2::track::PID::getMass(dautherPos); - float massNeg = o2::track::PID::getMass(dautherNeg); - // Axes specifications AxisSpec posZaxis = {400, -20., 20., "V_{z} (cm)"}; - AxisSpec dcaXYaxis = {120, -3.0, 3.0, "DCA_{xy} (cm)"}; - AxisSpec dcaZaxis = {120, -3.0, 3.0, "DCA_{z} (cm)"}; + AxisSpec dcaXYaxis = {800, -2.0, 2.0, "DCA_{xy} (cm)"}; + AxisSpec dcaZaxis = {800, -2.0, 2.0, "DCA_{z} (cm)"}; AxisSpec etaQAaxis = {1000, -1.0, 1.0, "#eta"}; AxisSpec tpcNClsFoundQAaxis = {110, 50., 160., "tpcNClsFound"}; @@ -109,45 +110,51 @@ struct PhianalysisTHnSparse { Service pdg; int n = 0; + float massPos = o2::track::PID::getMass(3); + float massNeg = o2::track::PID::getMass(3); double* pointPair = nullptr; double* pointSys = nullptr; - TLorentzVector d1, d2, mother; - bool produceQA, dataQA, MCTruthQA, t1, t2 = false; - int id; - float etapair = 0.8; - float tpcnSigmaPos = 3.0f; - float tpcnSigmaNeg = 3.0f; + ROOT::Math::PxPyPzMVector d1, d2, mother; + bool produceQA, dataQA, MCTruthQA, globalTrack, tpcPidOnly = false; + float tpcnSigmaPos = 100.0f; + float tpcnSigmaNeg = 100.0f; + float combinedNSigma = 100.0f; + float ptTOFThreshold = 0.5f; int tpcNClsFound = 70; + int dauSize = 2; rsn::MixingType mixingType = rsn::MixingType::none; + float vzCut = 10.0f; Filter triggerFilter = (o2::aod::evsel::sel8 == true); - Filter vtxFilter = (nabs(o2::aod::collision::posZ) < static_cast(cut.vZ)); - - Filter ptFilter = nabs(aod::track::pt) > static_cast(cut.pt); - Filter etaFilter = nabs(aod::track::eta) < static_cast(cut.etatrack); - Filter dcaFilter = (nabs(o2::aod::track::dcaXY) < static_cast(cut.dcaXY)) && (nabs(o2::aod::track::dcaZ) < static_cast(cut.dcaZ)); + Filter vtxFilter = (nabs(o2::aod::collision::posZ) < vzCut); using EventCandidates = soa::Filtered>; using EventCandidate = EventCandidates::iterator; - using TrackCandidates = soa::Filtered>; + using TrackCandidates = soa::Join; using EventCandidatesMC = soa::Filtered>; - using TrackCandidatesMC = soa::Filtered>; + using TrackCandidatesMC = soa::Join; using EventCandidatesMCGen = soa::Join; - using LabeledTracks = soa::Join; Preslice perCollision = aod::track::collisionId; - Partition positive = (aod::track::signed1Pt > 0.0f) && (nabs(o2::aod::pidtpc::tpcNSigmaKa) < std::abs(static_cast(cut.tpcnSigmaPos))); - Partition negative = (aod::track::signed1Pt < 0.0f) && (nabs(o2::aod::pidtpc::tpcNSigmaKa) < std::abs(static_cast(cut.tpcnSigmaNeg))); + Partition positive = (aod::track::signed1Pt > 0.0f); + Partition negative = (aod::track::signed1Pt < 0.0f); - Partition positiveMC = (aod::track::signed1Pt > 0.0f) && (nabs(o2::aod::pidtpc::tpcNSigmaKa) < std::abs(static_cast(cut.tpcnSigmaPos))); - Partition negativeMC = (aod::track::signed1Pt < 0.0f) && (nabs(o2::aod::pidtpc::tpcNSigmaKa) < std::abs(static_cast(cut.tpcnSigmaNeg))); + Partition positiveMC = (aod::track::signed1Pt > 0.0f); + Partition negativeMC = (aod::track::signed1Pt < 0.0f); void init(o2::framework::InitContext&) { + // defined in DataFormats/Reconstruction/include/ReconstructionDataFormats/PID.h + massPos = o2::track::PID::getMass(static_cast(daughterPos)); + massNeg = o2::track::PID::getMass(static_cast(daughterNeg)); + LOGF(info, "Initializing particle masses: "); + LOGF(info, " Positive: %d, mass: %f", static_cast(daughterPos), massPos); + LOGF(info, " Negative: %d, mass: %f", static_cast(daughterNeg), massNeg); + // Sparse axes AxisSpec invAxis = {invaxis, "Inv. mass (GeV/c^{2})", "im"}; AxisSpec ptAxis = {ptaxis, "p_{T} (GeV/c)", "pt"}; @@ -172,85 +179,113 @@ struct PhianalysisTHnSparse { produceQA = static_cast(produce.produceQA); mixingType = rsn::mixingTypeName(static_cast(produce.eventMixing)); - etapair = static_cast(cut.etapair); tpcnSigmaPos = static_cast(cut.tpcnSigmaPos); tpcnSigmaNeg = static_cast(cut.tpcnSigmaNeg); tpcNClsFound = static_cast(cut.tpcNClsFound); + globalTrack = static_cast(cut.globalTrack); + combinedNSigma = static_cast(cut.combinedNSigma); + tpcPidOnly = static_cast(cut.tpcPidOnly); + ptTOFThreshold = static_cast(cut.ptTOFThreshold); pointPair = new double[static_cast(o2::analysis::rsn::PairAxisType::unknown)]; pointSys = new double[static_cast(o2::analysis::rsn::SystematicsAxisType::unknown)]; rsnOutput = new o2::analysis::rsn::OutputSparse(); rsnOutput->init(sparseAxes, allAxes, sysAxes, allAxesSys, static_cast(produce.produceTrue), mixingType, static_cast(produce.produceLikesign), ®istry); + // Print summary of configuration + LOGF(info, "=== PhianalysisTHnSparse configuration summary ==="); + LOGF(info, "produceQA: %s", produceQA ? "true" : "false"); + LOGF(info, "produceTrue: %s", static_cast(produce.produceTrue) ? "true" : "false"); + LOGF(info, "produceLikesign: %s", static_cast(produce.produceLikesign) ? "true" : "false"); + LOGF(info, "eventMixing: %s", static_cast(produce.eventMixing).c_str()); + LOGF(info, "daughterPos: %d (PDG: %d)", static_cast(daughterPos), static_cast(daughterPosPDG)); + LOGF(info, "daughterNeg: %d (PDG: %d)", static_cast(daughterNeg), static_cast(daughterNegPDG)); + LOGF(info, "motherPDG: %d", static_cast(motherPDG)); + LOGF(info, "tpcnSigmaPos: %.2f", tpcnSigmaPos); + LOGF(info, "tpcnSigmaNeg: %.2f", tpcnSigmaNeg); + LOGF(info, "tpcPidOnly: %s", tpcPidOnly ? "true" : "false"); + LOGF(info, "combinedNSigma: %.2f", combinedNSigma); + LOGF(info, "ptTOFThreshold: %.2f", ptTOFThreshold); + LOGF(info, "rapidity: %.2f", static_cast(cut.rapidity)); + LOGF(info, "etatrack: %.2f", static_cast(cut.etatrack)); + LOGF(info, "etapair: %.2f", static_cast(cut.etapair)); + LOGF(info, "pt (min): %.2f", static_cast(cut.pt)); + LOGF(info, "dcaXY: %.2f", static_cast(cut.dcaXY)); + LOGF(info, "dcaZ: %.2f", static_cast(cut.dcaZ)); + LOGF(info, "globalTrack: %s", globalTrack ? "true" : "false"); + LOGF(info, "tpcNClsFound: %d", tpcNClsFound); + LOGF(info, "mixingType: %d", static_cast(mixingType)); + LOGF(info, "numberofMixedEvents: %d", static_cast(numberofMixedEvents)); + LOGF(info, "sparseAxes: "); + for (const auto& axis : static_cast>(sparseAxes)) { + LOGF(info, " %s", axis.c_str()); + } + LOGF(info, "sysAxes: "); + for (const auto& axis : static_cast>(sysAxes)) { + LOGF(info, " %s", axis.c_str()); + } + LOGF(info, "==============================================="); + if (produceQA) { // Event QA - registry.add("QAEvent/hSelection", "Event selection statistics", kTH1F, {{1, 0.0f, 1.0f}}); + registry.add("QAEvent/hSelection", "Event selection statistics", kTH1D, {{3, 0.0f, 3.0f}}); auto hEvent = registry.get(HIST("QAEvent/hSelection")); hEvent->GetXaxis()->SetBinLabel(1, "Full event statistics"); + hEvent->GetXaxis()->SetBinLabel(2, "Events with at least one #phi candidate"); + hEvent->GetXaxis()->SetBinLabel(3, "#phi candidates"); hEvent->SetMinimum(0.1); registry.add("QAEvent/hVtxZ", "Vertex position along the z-axis", kTH1F, {posZaxis}); - registry.add("QAEvent/hCent", "Distribution of multiplicity percentile", kTH1F, {{100, 0., 100.}}); + registry.add("QAEvent/hCent", "Distribution of multiplicity percentile", kTH1F, {{101, 0., 101.}}); registry.add("QAEvent/hMult", "Multiplicity (amplitude of non-zero channels in the FV0A + FV0C) ", kTH1F, {{300, 0., 30000.}}); - registry.add("QAEvent/h2Size", "Number of positive vs. negative Kaons per collision", kTH2F, {{30, 0., 30.}, {30, 0., 30.}}); // Track QA - registry.add("QATrack/hSelection", "Tracks combinations statistics", kTH1F, {{11, 0.0f, 11.0f}}); + registry.add("QATrack/hSelection", "Tracks statistics", kTH1D, {{10, 0.0f, 10.0f}}); auto hTrack = registry.get(HIST("QATrack/hSelection")); - hTrack->GetXaxis()->SetBinLabel(1, "all K^{+} K^{-} combinations"); - hTrack->GetXaxis()->SetBinLabel(2, "all K^{+}"); - hTrack->GetXaxis()->SetBinLabel(3, "all K^{-}"); - hTrack->GetXaxis()->SetBinLabel(4, "K^{+} tpcNClsFound"); - hTrack->GetXaxis()->SetBinLabel(5, "K^{-} tpcNClsFound"); - hTrack->GetXaxis()->SetBinLabel(6, "K^{+} isPrimaryTrack"); - hTrack->GetXaxis()->SetBinLabel(7, "K^{-} isPrimaryTrack"); - hTrack->GetXaxis()->SetBinLabel(8, "K^{+} isPVContributor"); - hTrack->GetXaxis()->SetBinLabel(9, "K^{-} isPVContributor"); - hTrack->GetXaxis()->SetBinLabel(10, "selected combinations"); - hTrack->GetXaxis()->SetBinLabel(11, "selected pairs (eta cut)"); + hTrack->GetXaxis()->SetBinLabel(1, "all tracks"); + hTrack->GetXaxis()->SetBinLabel(2, "passed pT cut"); + hTrack->GetXaxis()->SetBinLabel(3, "passed eta cut"); + hTrack->GetXaxis()->SetBinLabel(4, "passed DCA cut"); + hTrack->GetXaxis()->SetBinLabel(5, "passed PID cut"); + hTrack->GetXaxis()->SetBinLabel(6, "passed rapidity cut"); + hTrack->GetXaxis()->SetBinLabel(7, "passed tpcNClsFound cut"); + hTrack->GetXaxis()->SetBinLabel(8, "passed isPrimaryTrack cut"); + hTrack->GetXaxis()->SetBinLabel(9, "passed isPVContributor cut"); + hTrack->GetXaxis()->SetBinLabel(10, "passed all cuts"); hTrack->SetMinimum(0.1); - // Track1 - registry.add("QATrack/bs/hTrack1pt", "K^{+} p_{T} before selection", kTH1F, {ptAxis}); - registry.add("QATrack/bs/hTrack1DCAxy", "K^{+} DCA_{xy} before selection", kTH1F, {dcaXYaxis}); - registry.add("QATrack/bs/hTrack1DCAz", "K^{+} DCA_{z} before selection", kTH1F, {dcaZaxis}); - registry.add("QATrack/bs/hTrack1eta", "K^{+} #eta before selection", kTH1F, {{etaQAaxis}}); - registry.add("QATrack/bs/hTrack1tpcNClsFound", "K^{+} tpcNClsFound before selection", kTH1F, {{tpcNClsFoundQAaxis}}); - - registry.add("QATrack/as/hTrack1pt", "K^{+} p_{T} after selection", kTH1F, {ptAxis}); - registry.add("QATrack/as/hTrack1DCAxy", "K^{+} DCA_{xy} after selection", kTH1F, {dcaXYaxis}); - registry.add("QATrack/as/hTrack1DCAz", "K^{+} DCA_{z} after selection", kTH1F, {dcaZaxis}); - registry.add("QATrack/as/hTrack1eta", "K^{+} #eta after selection", kTH1F, {{etaQAaxis}}); - registry.add("QATrack/as/hTrack1tpcNClsFound", "K^{+} tpcNClsFound after selection", kTH1F, {{tpcNClsFoundQAaxis}}); - - // Track2 - registry.add("QATrack/bs/hTrack2pt", "K^{-} p_{T} before selection", kTH1F, {ptAxis}); - registry.add("QATrack/bs/hTrack2DCAxy", "K^{-} DCA_{xy} before selection", kTH1F, {dcaXYaxis}); - registry.add("QATrack/bs/hTrack2DCAz", "K^{-} DCA_{z} before selection", kTH1F, {dcaZaxis}); - registry.add("QATrack/bs/hTrack2eta", "K^{-} #eta before selection", kTH1F, {{etaQAaxis}}); - registry.add("QATrack/bs/hTrack2tpcNClsFound", "K^{-} tpcNClsFound before selection", kTH1F, {{tpcNClsFoundQAaxis}}); - - registry.add("QATrack/as/hTrack2pt", "K^{-} p_{T} after selection", kTH1F, {ptAxis}); - registry.add("QATrack/as/hTrack2DCAxy", "K^{-} DCA_{xy} after selection", kTH1F, {dcaXYaxis}); - registry.add("QATrack/as/hTrack2DCAz", "K^{-} DCA_{z} after selection", kTH1F, {dcaZaxis}); - registry.add("QATrack/as/hTrack2eta", "K^{-} #eta after selection", kTH1F, {{etaQAaxis}}); - registry.add("QATrack/as/hTrack2tpcNClsFound", "K^{-} tpcNClsFound after selection", kTH1F, {{tpcNClsFoundQAaxis}}); + registry.add("QATrack/hRapidity", "Rapidity distribution of K^{+} and K^{-}", kTH1F, {{200, -1, 1}}); + registry.add("QATrack/hEta", "Pseudorapidity distribution of K^{+} and K^{-}", kTH1F, {{200, -1, 1}}); + registry.add("QATrack/hTPCNClsFound", "Distribution of TPC NClsFound of K^{+} and K^{-}", kTH1F, {tpcNClsFoundQAaxis}); + registry.add("QATrack/hDCAxy", "Distribution of DCA_{xy} of K^{+} and K^{-}", kTH1F, {dcaXYaxis}); + registry.add("QATrack/hDCAz", "Distribution of DCA_{z} of K^{+} and K^{-}", kTH1F, {dcaZaxis}); + registry.add("QATrack/hPt", "Distribution of p_{T} of K^{+} and K^{-}", kTH1F, {ptaxis}); // TPC PID - registry.add("QATrack/TPCPID/h2TracknSigma", "", kTH2F, {{120, -6, 6}, {120, -6, 6}}); - auto h2TracknSigma = registry.get(HIST("QATrack/TPCPID/h2TracknSigma")); - h2TracknSigma->GetXaxis()->SetTitle("n#sigma_{TPC} K^{+}"); - h2TracknSigma->GetYaxis()->SetTitle("n#sigma_{TPC} K^{-}"); + registry.add("QATrack/hTPCnSigma", "Distribution of TPC nSigma of K^{+} and K^{-}", kTH1F, {{200, -10, 10}}); - registry.add("QATrack/TPCPID/h2nTrack1SigmaPt", "", kTH2F, {{100, 0, 10}, {120, -6, 6}}); - auto h2nTrack1SigmaPt = registry.get(HIST("QATrack/TPCPID/h2nTrack1SigmaPt")); - h2nTrack1SigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); - h2nTrack1SigmaPt->GetYaxis()->SetTitle("n#sigma_{TPC} K^{+}"); + registry.add("QATrack/h2TPCnSigma", "", kTH2F, {{200, -10, 10}, {200, -10, 10}}); + auto h2TPCnSigma = registry.get(HIST("QATrack/h2TPCnSigma")); + h2TPCnSigma->GetXaxis()->SetTitle("n#sigma_{TPC} K^{+}"); + h2TPCnSigma->GetYaxis()->SetTitle("n#sigma_{TPC} K^{-}"); - registry.add("QATrack/TPCPID/h2nTrack2SigmaPt", "", kTH2F, {{100, 0, 10}, {120, -6, 6}}); - auto h2nTrack2SigmaPt = registry.get(HIST("QATrack/TPCPID/h2nTrack2SigmaPt")); - h2nTrack2SigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); - h2nTrack2SigmaPt->GetYaxis()->SetTitle("n#sigma_{TPC} K^{-}"); + registry.add("QATrack/h2TPCnSigmaPt", "", kTH2F, {ptaxis, {200, -10, 10}}); + auto h2TPCnSigmaPt = registry.get(HIST("QATrack/h2TPCnSigmaPt")); + h2TPCnSigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + h2TPCnSigmaPt->GetYaxis()->SetTitle("n#sigma_{TPC} K^{#pm}"); + + // TOF PID + registry.add("QATrack/hTOFnSigma", "Distribution of TOF nSigma of K^{+} and K^{-}", kTH1F, {{200, -10, 10}}); + + registry.add("QATrack/h2TOFnSigma", "", kTH2F, {{200, -10, 10}, {200, -10, 10}}); + auto h2TOFnSigma = registry.get(HIST("QATrack/h2TOFnSigma")); + h2TOFnSigma->GetXaxis()->SetTitle("n#sigma_{TOF} K^{+}"); + h2TOFnSigma->GetYaxis()->SetTitle("n#sigma_{TOF} K^{-}"); + + registry.add("QATrack/h2TOFnSigmaPt", "", kTH2F, {ptaxis, {200, -10, 10}}); + auto h2TOFnSigmaPt = registry.get(HIST("QATrack/h2TOFnSigmaPt")); + h2TOFnSigmaPt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + h2TOFnSigmaPt->GetYaxis()->SetTitle("n#sigma_{TOF} K^{#pm}"); // MC Truth if (static_cast(produce.produceTrue)) { @@ -259,27 +294,6 @@ struct PhianalysisTHnSparse { hMCEventTruth->GetXaxis()->SetBinLabel(1, "Full MC Truth event statistics"); hMCEventTruth->SetMinimum(0.1); - registry.add("QAMC/Truth/hMCTrack", "MC Truth Track statistics", kTH1F, {{17, 0.0f, 17.0f}}); - auto hMCTrackTruth = registry.get(HIST("QAMC/Truth/hMCTrack")); - hMCTrackTruth->GetXaxis()->SetBinLabel(1, "all K^{+} K^{-} combinations"); - hMCTrackTruth->GetXaxis()->SetBinLabel(2, "all K^{+}"); - hMCTrackTruth->GetXaxis()->SetBinLabel(3, "all K^{-}"); - hMCTrackTruth->GetXaxis()->SetBinLabel(4, "K^{+} tpcNClsFound"); - hMCTrackTruth->GetXaxis()->SetBinLabel(5, "K^{-} tpcNClsFound"); - hMCTrackTruth->GetXaxis()->SetBinLabel(6, "K^{+} isPrimaryTrack"); - hMCTrackTruth->GetXaxis()->SetBinLabel(7, "K^{-} isPrimaryTrack"); - hMCTrackTruth->GetXaxis()->SetBinLabel(8, "K^{+} isPVContributor"); - hMCTrackTruth->GetXaxis()->SetBinLabel(9, "K^{-} isPVContributor"); - hMCTrackTruth->GetXaxis()->SetBinLabel(10, "selected combinations"); - hMCTrackTruth->GetXaxis()->SetBinLabel(11, "MCtrack PDG = 321"); - hMCTrackTruth->GetXaxis()->SetBinLabel(12, "all mothers"); - hMCTrackTruth->GetXaxis()->SetBinLabel(13, "equal mother PDGs"); - hMCTrackTruth->GetXaxis()->SetBinLabel(14, "equal mother IDs"); - hMCTrackTruth->GetXaxis()->SetBinLabel(15, "mother rapidity cut"); - hMCTrackTruth->GetXaxis()->SetBinLabel(16, "mother PDG = 333"); - hMCTrackTruth->GetXaxis()->SetBinLabel(17, "selected pairs (eta cut)"); - hMCTrackTruth->SetMinimum(0.1); - registry.add("QAMC/hInvMassTrueFalse", "", kTH1F, {invAxis}); // not written events in True distribution due to repetition of mothers?? // MC Gen @@ -289,33 +303,14 @@ struct PhianalysisTHnSparse { hMCEventGen->GetXaxis()->SetBinLabel(2, "McCollision V_{z} cut"); hMCEventGen->GetXaxis()->SetBinLabel(3, "collisions"); hMCEventGen->SetMinimum(0.1); - - registry.add("QAMC/Gen/hMCTrack", "MC Gen Track statistics", kTH1D, {{7, 0.0f, 7.0f}}); - auto hMCTrackGen = registry.get(HIST("QAMC/Gen/hMCTrack")); - hMCTrackGen->GetXaxis()->SetBinLabel(1, "all mcParticles"); - hMCTrackGen->GetXaxis()->SetBinLabel(2, "rapidity cut"); - hMCTrackGen->GetXaxis()->SetBinLabel(3, "particle PDG = 333"); - hMCTrackGen->GetXaxis()->SetBinLabel(4, "has 2 dauthers"); - hMCTrackGen->GetXaxis()->SetBinLabel(5, "all dauthers"); - hMCTrackGen->GetXaxis()->SetBinLabel(6, "isPhysicalPrimary"); - hMCTrackGen->GetXaxis()->SetBinLabel(7, "selected pairs"); - hMCTrackGen->SetMinimum(0.1); } // Mixing QA if (mixingType != rsn::MixingType::none) { - registry.add("QAMixing/hSelection", "Event mixing selection statistics", kTH1F, {{1, 0.0f, 1.0f}}); + registry.add("QAMixing/hSelection", "Event mixing selection statistics", kTH1D, {{1, 0.0f, 1.0f}}); auto hEM = registry.get(HIST("QAMixing/hSelection")); hEM->GetXaxis()->SetBinLabel(1, "Full event mixing statistics"); hEM->SetMinimum(0.1); - registry.add("QAMixing/hTrackSelection", "Event mixing tracks combinations statistics", kTH1F, {{4, 0.0f, 4.0f}}); - auto hEMTrack = registry.get(HIST("QAMixing/hTrackSelection")); - hEMTrack->GetXaxis()->SetBinLabel(1, "all K^{+} K^{-} combinations"); - hEMTrack->GetXaxis()->SetBinLabel(2, "all K^{+}"); - hEMTrack->GetXaxis()->SetBinLabel(3, "all K^{-}"); - hEMTrack->GetXaxis()->SetBinLabel(4, "selected pairs (eta cut)"); - hEMTrack->SetMinimum(0.1); - registry.add("QAMixing/h2mu1_mu2", "Event Mixing Multiplicity", kTH2F, {axisMultiplicityMixing, axisMultiplicityMixing}); auto h2EMmu = registry.get(HIST("QAMixing/h2mu1_mu2")); h2EMmu->GetXaxis()->SetTitle("1.Event multiplicity"); @@ -338,88 +333,94 @@ struct PhianalysisTHnSparse { } template - bool selectedTrack(const T& track) + bool selectedTrack(const T& track, bool isPositive) { - if (produceQA) { - if (t1) { - if (dataQA) - registry.fill(HIST("QATrack/hSelection"), 1.5); - if (MCTruthQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 1.5); - } - if (t2) { - if (dataQA) - registry.fill(HIST("QATrack/hSelection"), 2.5); - if (MCTruthQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 2.5); - } + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 0.5); // all tracks + + // Apply pT cut + if (track.pt() < static_cast(cut.pt)) + return false; + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 1.5); + + // Apply eta cut + if (std::abs(track.eta()) >= static_cast(cut.etatrack)) + return false; + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 2.5); + + // Apply DCA cuts + if (std::abs(track.dcaXY()) >= static_cast(cut.dcaXY) || + std::abs(track.dcaZ()) >= static_cast(cut.dcaZ)) + return false; + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 3.5); + + // PID selection: TPC-only for pt < threshold value, TPC+TOF for pt >= threshold value and have TOF, else TPC-only + float nSigmaCut = isPositive ? tpcnSigmaPos : tpcnSigmaNeg; + if (track.pt() < ptTOFThreshold || !track.hasTOF() || tpcPidOnly) { + if (std::abs(track.tpcNSigmaKa()) >= nSigmaCut) + return false; + } else { + if (std::sqrt(track.tpcNSigmaKa() * track.tpcNSigmaKa() + track.tofNSigmaKa() * track.tofNSigmaKa()) >= combinedNSigma) + return false; } + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 4.5); - if (track.tpcNClsFound() < tpcNClsFound) + // Apply rapidity cut + if (std::abs(track.rapidity(isPositive ? massPos : massNeg)) > static_cast(cut.rapidity)) { return false; - if (produceQA) { - if (t1) { - if (dataQA) - registry.fill(HIST("QATrack/hSelection"), 3.5); - if (MCTruthQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 3.5); - } - if (t2) { - if (dataQA) - registry.fill(HIST("QATrack/hSelection"), 4.5); - if (MCTruthQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 4.5); - } } - if (!track.isPrimaryTrack()) + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 5.5); + + // Apply tpcNClsFound cut + if (track.tpcNClsFound() < tpcNClsFound) return false; - if (produceQA) { - if (t1) { - if (dataQA) - registry.fill(HIST("QATrack/hSelection"), 5.5); - if (MCTruthQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 5.5); - } - if (t2) { - if (dataQA) - registry.fill(HIST("QATrack/hSelection"), 6.5); - if (MCTruthQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 6.5); - } + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 6.5); + + if (globalTrack) { + // Apply Global track cuts + if (!track.isGlobalTrack()) + return false; + } else { + // Apply Primary track cuts + if (!track.isPrimaryTrack()) + return false; } + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 7.5); + + // Apply PV Contributor cuts if (!track.isPVContributor()) return false; - if (produceQA) { - if (t1) { - if (dataQA) - registry.fill(HIST("QATrack/hSelection"), 7.5); - if (MCTruthQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 7.5); - } - if (t2) { - if (dataQA) - registry.fill(HIST("QATrack/hSelection"), 8.5); - if (MCTruthQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 8.5); - } - } + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 8.5); + + if (produceQA && dataQA) + registry.fill(HIST("QATrack/hSelection"), 9.5); + return true; } template - bool selectedPair(TLorentzVector& mother, const T& track1, const T& track2) + bool selectedPair(ROOT::Math::PxPyPzMVector& mother, const T& track1, const T& track2) { - d1.SetXYZM(track1.px(), track1.py(), track1.pz(), massPos); - d2.SetXYZM(track2.px(), track2.py(), track2.pz(), massNeg); + d1 = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massPos); + d2 = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massNeg); mother = d1 + d2; - if (std::abs(mother.Eta()) > etapair) + if (std::abs(mother.Eta()) > static_cast(cut.etapair)) return false; + return true; } template float getMultiplicity(const T& collision) { - float multiplicity = collision.multFT0C() + collision.multFT0A(); + float multiplicity = collision.multFT0M(); return multiplicity; } template @@ -449,81 +450,74 @@ struct PhianalysisTHnSparse { void processData(EventCandidate const& collision, TrackCandidates const& /*tracks*/) { - auto posDauthers = positive->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); - auto negDauthers = negative->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); - + auto posDaughters = positive->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto negDaughters = negative->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + n = 0; if (produceQA) { registry.fill(HIST("QAEvent/hSelection"), 0.5); - registry.fill(HIST("QAEvent/h2Size"), posDauthers.size(), negDauthers.size()); registry.fill(HIST("QAEvent/hVtxZ"), collision.posZ()); registry.fill(HIST("QAEvent/hMult"), getMultiplicity(collision)); registry.fill(HIST("QAEvent/hCent"), getCentrality(collision)); + + dataQA = true; + for (const auto& track : posDaughters) { + selectedTrack(track, true); + } + for (const auto& track : negDaughters) { + selectedTrack(track, false); + } + dataQA = false; } if (static_cast(verbose.verboselevel) > 0 && static_cast(verbose.refresh) > 0 && collision.globalIndex() % static_cast(verbose.refresh) == static_cast(verbose.refreshIndex)) - LOGF(info, "%d pos=%lld neg=%lld, Z vertex position: %f [cm]", collision.globalIndex(), posDauthers.size(), negDauthers.size(), collision.posZ()); + LOGF(info, "%d pos=%lld neg=%lld, Z vertex position: %f [cm]", collision.globalIndex(), posDaughters.size(), negDaughters.size(), collision.posZ()); - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthers, negDauthers))) { - if (produceQA) { - registry.fill(HIST("QATrack/hSelection"), 0.5); - - registry.fill(HIST("QATrack/bs/hTrack1pt"), track1.pt()); - registry.fill(HIST("QATrack/bs/hTrack1DCAxy"), track1.dcaXY()); - registry.fill(HIST("QATrack/bs/hTrack1DCAz"), track1.dcaZ()); - registry.fill(HIST("QATrack/bs/hTrack1eta"), track1.eta()); - registry.fill(HIST("QATrack/bs/hTrack1tpcNClsFound"), track1.tpcNClsFound()); - - registry.fill(HIST("QATrack/bs/hTrack2pt"), track2.pt()); - registry.fill(HIST("QATrack/bs/hTrack2DCAxy"), track2.dcaXY()); - registry.fill(HIST("QATrack/bs/hTrack2DCAz"), track2.dcaZ()); - registry.fill(HIST("QATrack/bs/hTrack2eta"), track2.eta()); - registry.fill(HIST("QATrack/bs/hTrack2tpcNClsFound"), track2.tpcNClsFound()); - } + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughters, negDaughters))) { - dataQA = true; - t1 = true; - if (!selectedTrack(track1)) + if (!selectedTrack(track1, true)) // track1 is positive continue; - t1 = false; - t2 = true; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, false)) // track2 is negative continue; - t2 = false; - dataQA = false; - - if (produceQA) { - registry.fill(HIST("QATrack/hSelection"), 9.5); - - registry.fill(HIST("QATrack/as/hTrack1pt"), track1.pt()); - registry.fill(HIST("QATrack/as/hTrack1DCAxy"), track1.dcaXY()); - registry.fill(HIST("QATrack/as/hTrack1DCAz"), track1.dcaZ()); - registry.fill(HIST("QATrack/as/hTrack1eta"), track1.eta()); - registry.fill(HIST("QATrack/as/hTrack1tpcNClsFound"), track1.tpcNClsFound()); - - registry.fill(HIST("QATrack/as/hTrack2pt"), track2.pt()); - registry.fill(HIST("QATrack/as/hTrack2DCAxy"), track2.dcaXY()); - registry.fill(HIST("QATrack/as/hTrack2DCAz"), track2.dcaZ()); - registry.fill(HIST("QATrack/as/hTrack2eta"), track2.eta()); - registry.fill(HIST("QATrack/as/hTrack2tpcNClsFound"), track2.tpcNClsFound()); - } if (!selectedPair(mother, track1, track2)) continue; if (produceQA) { - registry.fill(HIST("QATrack/hSelection"), 10.5); - - registry.fill(HIST("QATrack/TPCPID/h2TracknSigma"), track1.tpcNSigmaKa(), track2.tpcNSigmaKa()); - registry.fill(HIST("QATrack/TPCPID/h2nTrack1SigmaPt"), track1.pt(), track1.tpcNSigmaKa()); - registry.fill(HIST("QATrack/TPCPID/h2nTrack2SigmaPt"), track2.pt(), track2.tpcNSigmaKa()); + registry.fill(HIST("QATrack/h2TPCnSigma"), track1.tpcNSigmaKa(), track2.tpcNSigmaKa()); + registry.fill(HIST("QATrack/h2TPCnSigmaPt"), track1.pt(), track1.tpcNSigmaKa()); + registry.fill(HIST("QATrack/h2TPCnSigmaPt"), track2.pt(), track2.tpcNSigmaKa()); + + registry.fill(HIST("QATrack/h2TOFnSigma"), track1.tofNSigmaKa(), track2.tofNSigmaKa()); + registry.fill(HIST("QATrack/h2TOFnSigmaPt"), track1.pt(), track1.tofNSigmaKa()); + registry.fill(HIST("QATrack/h2TOFnSigmaPt"), track2.pt(), track2.tofNSigmaKa()); + + registry.fill(HIST("QATrack/hTPCnSigma"), track1.tpcNSigmaKa()); + registry.fill(HIST("QATrack/hTPCnSigma"), track2.tpcNSigmaKa()); + if (track1.hasTOF()) + registry.fill(HIST("QATrack/hTOFnSigma"), track1.tofNSigmaKa()); + if (track2.hasTOF()) + registry.fill(HIST("QATrack/hTOFnSigma"), track2.tofNSigmaKa()); + + registry.fill(HIST("QATrack/hEta"), track1.eta()); + registry.fill(HIST("QATrack/hEta"), track2.eta()); + registry.fill(HIST("QATrack/hPt"), track1.pt()); + registry.fill(HIST("QATrack/hPt"), track2.pt()); + registry.fill(HIST("QATrack/hDCAxy"), track1.dcaXY()); + registry.fill(HIST("QATrack/hDCAxy"), track2.dcaXY()); + registry.fill(HIST("QATrack/hDCAz"), track1.dcaZ()); + registry.fill(HIST("QATrack/hDCAz"), track2.dcaZ()); + registry.fill(HIST("QATrack/hTPCNClsFound"), track1.tpcNClsFound()); + registry.fill(HIST("QATrack/hTPCNClsFound"), track2.tpcNClsFound()); + registry.fill(HIST("QATrack/hRapidity"), track1.rapidity(massPos)); + registry.fill(HIST("QATrack/hRapidity"), track2.rapidity(massNeg)); } - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(collision), getCentrality(collision), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), collision.posZ(), @@ -531,28 +525,34 @@ struct PhianalysisTHnSparse { 0, 0); rsnOutput->fillUnlikepm(pointPair); + + if (produceQA) + registry.fill(HIST("QAEvent/hSelection"), 2.5); + if (n == 0) + registry.fill(HIST("QAEvent/hSelection"), 1.5); + n = n + 1; } if (static_cast(produce.produceLikesign)) { - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(posDauthers, posDauthers))) { - if (!selectedTrack(track1)) + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(posDaughters, posDaughters))) { + if (!selectedTrack(track1, true)) // both positive continue; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, true)) // both positive continue; if (!selectedPair(mother, track1, track2)) continue; if (static_cast(verbose.verboselevel) > 1) - LOGF(info, "Like-sign positive: d1=%ld , d2=%ld , mother=%f", track1.globalIndex(), track2.globalIndex(), mother.Mag()); + LOGF(info, "Like-sign positive: d1=%ld , d2=%ld , mother=%f", track1.globalIndex(), track2.globalIndex(), mother.M()); - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(collision), getCentrality(collision), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), collision.posZ(), @@ -563,24 +563,24 @@ struct PhianalysisTHnSparse { rsnOutput->fillLikepp(pointPair); } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(negDauthers, negDauthers))) { - if (!selectedTrack(track1)) + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(negDaughters, negDaughters))) { + if (!selectedTrack(track1, false)) // both negative continue; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, false)) // both negative continue; if (!selectedPair(mother, track1, track2)) continue; if (static_cast(verbose.verboselevel) > 1) - LOGF(info, "Like-sign negative: d1=%ld , d2=%ld , mother=%f", track1.globalIndex(), track2.globalIndex(), mother.Mag()); + LOGF(info, "Like-sign negative: d1=%ld , d2=%ld , mother=%f", track1.globalIndex(), track2.globalIndex(), mother.M()); - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(collision), getCentrality(collision), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), collision.posZ(), @@ -601,76 +601,57 @@ struct PhianalysisTHnSparse { registry.fill(HIST("QAMC/Truth/hMCEvent"), 0.5); - auto posDauthersMC = positiveMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); - auto negDauthersMC = negativeMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto posDaughtersMC = positiveMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto negDaughtersMC = negativeMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); if (!collision.has_mcCollision()) { - LOGF(warning, "No MC collision for this collision, skip..."); + if (static_cast(verbose.verboselevel) > 0) + LOGF(warning, "No MC collision for this collision, skip..."); return; } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersMC, negDauthersMC))) { - registry.fill(HIST("QAMC/Truth/hMCTrack"), 0.5); + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersMC, negDaughtersMC))) { if (!track1.has_mcParticle()) { - LOGF(warning, "No MC particle for track, skip..."); + if (static_cast(verbose.verboselevel) > 0) + LOGF(warning, "No MC particle for track, skip..."); continue; } if (!track2.has_mcParticle()) { - LOGF(warning, "No MC particle for track, skip..."); + if (static_cast(verbose.verboselevel) > 0) + LOGF(warning, "No MC particle for track, skip..."); continue; } - MCTruthQA = true; - t1 = true; - if (!selectedTrack(track1)) + if (!selectedTrack(track1, true)) // track1 is positive continue; - t1 = false; - t2 = true; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, false)) // track2 is negative continue; - t2 = false; - MCTruthQA = false; - - registry.fill(HIST("QAMC/Truth/hMCTrack"), 9.5); const auto mctrack1 = track1.mcParticle(); const auto mctrack2 = track2.mcParticle(); int track1PDG = std::abs(mctrack1.pdgCode()); int track2PDG = std::abs(mctrack2.pdgCode()); - if (!(track1PDG == dautherPosPDG && track2PDG == dautherNegPDG)) { + if (!(track1PDG == daughterPosPDG && track2PDG == daughterNegPDG)) { continue; } - if (produceQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 10.5); - n = 0; for (const auto& mothertrack1 : mctrack1.mothers_as()) { for (const auto& mothertrack2 : mctrack2.mothers_as()) { - if (produceQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 11.5); if (mothertrack1.pdgCode() != mothertrack2.pdgCode()) continue; - if (produceQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 12.5); if (mothertrack1.globalIndex() != mothertrack2.globalIndex()) continue; - if (produceQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 13.5); - if (std::abs(mothertrack1.y()) > static_cast(cut.y)) + if (std::abs(mothertrack1.y()) > static_cast(cut.rapidity)) continue; - if (produceQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 14.5); if (std::abs(mothertrack1.pdgCode()) != motherPDG) continue; - if (produceQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 15.5); if (static_cast(verbose.verboselevel) > 1) { LOGF(info, "True: %d, d1=%d (%ld), d2=%d (%ld), mother=%d (%ld)", n, mctrack1.pdgCode(), mctrack1.globalIndex(), mctrack2.pdgCode(), mctrack2.globalIndex(), mothertrack1.pdgCode(), mothertrack1.globalIndex()); @@ -679,21 +660,19 @@ struct PhianalysisTHnSparse { if (!selectedPair(mother, mctrack1, mctrack2)) continue; - if (produceQA) - registry.fill(HIST("QAMC/Truth/hMCTrack"), 16.5); if (n > 0) { if (produceQA) - registry.fill(HIST("QAMC/hInvMassTrueFalse"), mother.Mag()); + registry.fill(HIST("QAMC/hInvMassTrueFalse"), mother.M()); continue; } - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(collision), getCentrality(collision), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), collision.posZ(), @@ -711,76 +690,66 @@ struct PhianalysisTHnSparse { void processGen(aod::McCollision const& mcCollision, soa::SmallGroups const& collisions, LabeledTracks const& /*particles*/, aod::McParticles const& mcParticles) { - registry.fill(HIST("QAMC/Gen/hMCEvent"), 0.5); - if (std::abs(mcCollision.posZ()) > static_cast(cut.vZ)) + + if (produceQA) + registry.fill(HIST("QAMC/Gen/hMCEvent"), 0.5); + + if (std::abs(mcCollision.posZ()) > vzCut) return; - registry.fill(HIST("QAMC/Gen/hMCEvent"), 1.5); + + if (produceQA) + registry.fill(HIST("QAMC/Gen/hMCEvent"), 1.5); if (collisions.size() == 0) return; for (const auto& collision : collisions) { - registry.fill(HIST("QAMC/Gen/hMCEvent"), 2.5); + if (produceQA) + registry.fill(HIST("QAMC/Gen/hMCEvent"), 2.5); if (!collision.has_mcCollision()) { - LOGF(warning, "No McCollision for this collision, skip..."); + if (static_cast(verbose.verboselevel) > 0) + LOGF(warning, "No McCollision for this collision, skip..."); return; } - auto centralityGen = 0; - centralityGen = getCentrality(collision); - auto multiplicityGen = 0; - multiplicityGen = getMultiplicity(collision); + auto centralityGen = getCentrality(collision); + auto multiplicityGen = getMultiplicity(collision); for (const auto& particle : mcParticles) { - registry.fill(HIST("QAMC/Gen/hMCTrack"), 0.5); - if (std::abs(particle.y()) > static_cast(cut.y)) + if (std::abs(particle.y()) > static_cast(cut.rapidity)) continue; - registry.fill(HIST("QAMC/Gen/hMCTrack"), 1.5); - if (particle.pdgCode() == motherPDG) { - registry.fill(HIST("QAMC/Gen/hMCTrack"), 2.5); auto daughters = particle.daughters_as(); - if (daughters.size() != 2) + if (daughters.size() != dauSize) continue; - registry.fill(HIST("QAMC/Gen/hMCTrack"), 3.5); - auto daup = false; auto daun = false; for (const auto& dau : daughters) { - registry.fill(HIST("QAMC/Gen/hMCTrack"), 4.5); - - if (!dau.isPhysicalPrimary()) - continue; - - registry.fill(HIST("QAMC/Gen/hMCTrack"), 5.5); - - if (dau.pdgCode() == dautherPosPDG) { + if (dau.pdgCode() == daughterPosPDG) { daup = true; - d1.SetXYZM(dau.px(), dau.py(), dau.pz(), massPos); - } else if (dau.pdgCode() == -dautherNegPDG) { + d1 = ROOT::Math::PxPyPzMVector(dau.px(), dau.py(), dau.pz(), massPos); + } else if (dau.pdgCode() == -daughterNegPDG) { daun = true; - d2.SetXYZM(dau.px(), dau.py(), dau.pz(), massNeg); + d2 = ROOT::Math::PxPyPzMVector(dau.px(), dau.py(), dau.pz(), massNeg); } } - if (!daup && !daun) + if (!daup || !daun) continue; - registry.fill(HIST("QAMC/Gen/hMCTrack"), 6.5); - mother = d1 + d2; - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), multiplicityGen, centralityGen, - std::abs(static_cast(cut.tpcnSigmaPos) / 2.0), - std::abs(static_cast(cut.tpcnSigmaNeg) / 2.0), + 0, + 0, mother.Eta(), mother.Rapidity(), mcCollision.posZ(), @@ -793,7 +762,7 @@ struct PhianalysisTHnSparse { } } } - PROCESS_SWITCH(PhianalysisTHnSparse, processGen, "Process MC Mateched.", true); + PROCESS_SWITCH(PhianalysisTHnSparse, processGen, "Process MC Generated.", false); void processMixed(EventCandidates const& collisions, TrackCandidates const& tracks) { @@ -813,10 +782,10 @@ struct PhianalysisTHnSparse { if (produceQA) registry.fill(HIST("QAMixing/hSelection"), 0.5); - auto posDauthersc1 = positive->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); - auto posDauthersc2 = positive->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); - auto negDauthersc1 = negative->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); - auto negDauthersc2 = negative->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); + auto posDaughtersc1 = positive->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); + auto posDaughtersc2 = positive->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); + auto negDaughtersc1 = negative->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); + auto negDaughtersc2 = negative->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); if (produceQA) { registry.fill(HIST("QAMixing/h2mu1_mu2"), getMultiplicity(c1), getMultiplicity(c2)); @@ -824,30 +793,22 @@ struct PhianalysisTHnSparse { registry.fill(HIST("QAMixing/h2vz1_vz2"), c1.posZ(), c2.posZ()); } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersc1, negDauthersc2))) { - if (produceQA) - registry.fill(HIST("QAMixing/hTrackSelection"), 0.5); + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc1, negDaughtersc2))) { - if (!selectedTrack(track1)) + if (!selectedTrack(track1, true)) // track1 is positive continue; - if (produceQA) - registry.fill(HIST("QAMixing/hTrackSelection"), 1.5); - if (!selectedTrack(track2)) + if (!selectedTrack(track2, false)) // track2 is negative continue; - if (produceQA) - registry.fill(HIST("QAMixing/hTrackSelection"), 2.5); if (!selectedPair(mother, track1, track2)) continue; - if (produceQA) - registry.fill(HIST("QAMixing/hTrackSelection"), 3.5); - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(c1), getCentrality(c1), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), c1.posZ(), @@ -860,23 +821,22 @@ struct PhianalysisTHnSparse { if (static_cast(produce.produceLikesign)) { - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersc1, posDauthersc2))) { - - if (!selectedTrack(track1)) + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc1, posDaughtersc2))) { + if (!selectedTrack(track1, true)) // track1 is positive continue; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, true)) // track2 is positive continue; if (!selectedPair(mother, track1, track2)) continue; - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(c1), getCentrality(c1), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), c1.posZ(), @@ -887,22 +847,21 @@ struct PhianalysisTHnSparse { rsnOutput->fillMixingpp(pointPair); } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(negDauthersc1, negDauthersc2))) { - - if (!selectedTrack(track1)) + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(negDaughtersc1, negDaughtersc2))) { + if (!selectedTrack(track1, false)) continue; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, false)) continue; if (!selectedPair(mother, track1, track2)) continue; - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(c1), getCentrality(c1), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), c1.posZ(), @@ -914,23 +873,22 @@ struct PhianalysisTHnSparse { } } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersc2, negDauthersc1))) { - - if (!selectedTrack(track1)) + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc2, negDaughtersc1))) { + if (!selectedTrack(track1, true)) // track1 is positive continue; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, false)) // track2 is negative continue; if (!selectedPair(mother, track1, track2)) continue; - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(c1), getCentrality(c1), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), c1.posZ(), @@ -947,10 +905,10 @@ struct PhianalysisTHnSparse { if (produceQA) registry.fill(HIST("QAMixing/hSelection"), 0.5); - auto posDauthersc1 = positive->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); - auto posDauthersc2 = positive->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); - auto negDauthersc1 = negative->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); - auto negDauthersc2 = negative->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); + auto posDaughtersc1 = positive->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); + auto posDaughtersc2 = positive->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); + auto negDaughtersc1 = negative->sliceByCached(aod::track::collisionId, c1.globalIndex(), cache); + auto negDaughtersc2 = negative->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); if (produceQA) { registry.fill(HIST("QAMixing/h2mu1_mu2"), getMultiplicity(c1), getMultiplicity(c2)); @@ -958,30 +916,23 @@ struct PhianalysisTHnSparse { registry.fill(HIST("QAMixing/h2vz1_vz2"), c1.posZ(), c2.posZ()); } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersc1, negDauthersc2))) { - if (produceQA) - registry.fill(HIST("QAMixing/hTrackSelection"), 0.5); + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc1, negDaughtersc2))) { - if (!selectedTrack(track1)) + if (!selectedTrack(track1, true)) // track1 is positive continue; - if (produceQA) - registry.fill(HIST("QAMixing/hTrackSelection"), 1.5); - if (!selectedTrack(track2)) + + if (!selectedTrack(track2, false)) // track2 is negative continue; - if (produceQA) - registry.fill(HIST("QAMixing/hTrackSelection"), 2.5); if (!selectedPair(mother, track1, track2)) continue; - if (produceQA) - registry.fill(HIST("QAMixing/hTrackSelection"), 3.5); - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(c1), getCentrality(c1), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), c1.posZ(), @@ -994,23 +945,23 @@ struct PhianalysisTHnSparse { if (static_cast(produce.produceLikesign)) { - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersc1, posDauthersc2))) { + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc1, posDaughtersc2))) { - if (!selectedTrack(track1)) + if (!selectedTrack(track1, true)) // track1 is positive continue; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, true)) // track2 is positive continue; if (!selectedPair(mother, track1, track2)) continue; - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(c1), getCentrality(c1), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), c1.posZ(), @@ -1021,22 +972,22 @@ struct PhianalysisTHnSparse { rsnOutput->fillMixingpp(pointPair); } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(negDauthersc1, negDauthersc2))) { + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(negDaughtersc1, negDaughtersc2))) { - if (!selectedTrack(track1)) + if (!selectedTrack(track1, false)) continue; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, false)) continue; if (!selectedPair(mother, track1, track2)) continue; - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(c1), getCentrality(c1), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), c1.posZ(), @@ -1048,23 +999,23 @@ struct PhianalysisTHnSparse { } } - for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDauthersc2, negDauthersc1))) { + for (const auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(posDaughtersc2, negDaughtersc1))) { - if (!selectedTrack(track1)) + if (!selectedTrack(track1, true)) continue; - if (!selectedTrack(track2)) + if (!selectedTrack(track2, false)) continue; if (!selectedPair(mother, track1, track2)) continue; - pointPair = fillPointPair(mother.Mag(), + pointPair = fillPointPair(mother.M(), mother.Pt(), getMultiplicity(c1), getCentrality(c1), - (tpcnSigmaPos > 0) ? std::abs(track1.tpcNSigmaKa()) : track1.tpcNSigmaKa(), - (tpcnSigmaNeg > 0) ? std::abs(track2.tpcNSigmaKa()) : track2.tpcNSigmaKa(), + track1.tpcNSigmaKa(), + track2.tpcNSigmaKa(), mother.Eta(), mother.Rapidity(), c1.posZ(), @@ -1077,76 +1028,8 @@ struct PhianalysisTHnSparse { } } } - PROCESS_SWITCH(PhianalysisTHnSparse, processMixed, "Process Mixing Event.", true); - - void processGenOld(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles) - { - if (!static_cast(produce.produceTrue)) - return; - - registry.fill(HIST("QAMC/hMC"), 0.5); - - if (std::abs(mcCollision.posZ()) > static_cast(cut.vZ)) - return; - - registry.fill(HIST("QAMC/hMC"), 1.5); - - for (const auto& particle : mcParticles) { - registry.fill(HIST("QAMC/hMC"), 2.5); - if (std::abs(particle.y()) > static_cast(cut.y)) - continue; - - registry.fill(HIST("QAMC/hMC"), 3.5); - - if (particle.pdgCode() == motherPDG) { - auto daughters = particle.daughters_as(); - if (daughters.size() != 2) - continue; - - registry.fill(HIST("QAMC/hMC"), 4.5); - - auto daup = false; - auto daun = false; - - for (const auto& dau : daughters) { - if (!dau.isPhysicalPrimary()) - continue; - - if (dau.pdgCode() == dautherPosPDG) { - daup = true; - d1.SetXYZM(dau.px(), dau.py(), dau.pz(), massPos); - } else if (dau.pdgCode() == -dautherNegPDG) { - daun = true; - d2.SetXYZM(dau.px(), dau.py(), dau.pz(), massNeg); - } - } - if (!daup && !daun) - continue; - - registry.fill(HIST("QAMC/hMC"), 5.5); - - mother = d1 + d2; - - pointPair = fillPointPair(mother.Mag(), - mother.Pt(), - 0, - 0, - std::abs(static_cast(cut.tpcnSigmaPos) / 2.0), - std::abs(static_cast(cut.tpcnSigmaNeg) / 2.0), - mother.Eta(), - mother.Rapidity(), - mcCollision.posZ(), - 0, - 0, - 0); - - rsnOutput->fillUnlikegenOld(pointPair); - } - } - } - PROCESS_SWITCH(PhianalysisTHnSparse, processGenOld, "Process generated.", false); + PROCESS_SWITCH(PhianalysisTHnSparse, processMixed, "Process Mixing Event.", false); }; - WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{