diff --git a/PWGLF/TableProducer/Nuspex/particleCompositionCorrection.cxx b/PWGLF/TableProducer/Nuspex/particleCompositionCorrection.cxx index b2f5e201476..2538a80d930 100644 --- a/PWGLF/TableProducer/Nuspex/particleCompositionCorrection.cxx +++ b/PWGLF/TableProducer/Nuspex/particleCompositionCorrection.cxx @@ -38,6 +38,8 @@ #include #include +// #include "Framework/Logger.h" + using namespace o2; using namespace o2::framework; using namespace o2::ml; @@ -64,6 +66,7 @@ struct ParticleCompositionCorrection { Configurable skipAll{"skipAll", false, "run table producer in dummy mode, i.e. skip all computations and fill with 1"}; Configurable skipSec{"skipSec", false, "dont calculate weights for secondaries"}; + Configurable skipNeutralSec{"skipNeutralSec", true, "dont calculate weights for neutral secondaries"}; Configurable skipNonPhysicalPrim{"skipNonPhysicalPrim", true, "dont calculate weights for particles that are not (originating from) physical primaries; i.e. reject (decays of) pi0 etc."}; Configurable etaCut{"etaCut", 0.8f, "eta cut"}; Configurable ptMinCut{"ptMinCut", 0.15f, "pt min cut"}; @@ -118,6 +121,7 @@ void ParticleCompositionCorrection::init(InitContext const&) const int maxMult = 100; const int nBinsMult = maxMult + 1; const AxisSpec multAxis = {nBinsMult, -0.5, nBinsMult - 0.5, "#it{N}_{ch}", "mult"}; + const AxisSpec processAxis = {51, -0.5, 50.5, "process", "process"}; histos.add("multDist_INELg0", "", kTH1D, {multAxis}); histos.add("multDist_fid", "", kTH1D, {multAxis}); @@ -127,25 +131,15 @@ void ParticleCompositionCorrection::init(InitContext const&) histos.add("frac/data/proton", "", kTProfile, {ptAxis}); histos.add("frac/data/sigma", "", kTProfile, {ptAxis}); - histos.add("frac/data/kaon_mult", "", kTProfile2D, {multAxis, ptAxis}); - histos.add("frac/data/proton_mult", "", kTProfile2D, {multAxis, ptAxis}); - histos.add("frac/data/sigma_mult", "", kTProfile2D, {multAxis, ptAxis}); - histos.addClone("frac/data/", "frac/mc/"); histos.add("weight/kaon", "", kTProfile, {ptAxis}); histos.add("weight/proton", "", kTProfile, {ptAxis}); histos.add("weight/sigma", "", kTProfile, {ptAxis}); - histos.add("weight/kaon_mult", "", kTProfile2D, {multAxis, ptAxis}); - histos.add("weight/proton_mult", "", kTProfile2D, {multAxis, ptAxis}); - histos.add("weight/sigma_mult", "", kTProfile2D, {multAxis, ptAxis}); - + histos.add("weight/prim", "", kTProfile, {ptAxis}); histos.add("weight/secDec", "", kTProfile, {ptAxis}); histos.add("weight/secMat", "", kTProfile, {ptAxis}); - - histos.add("weight/secDec_mult", "", kTProfile2D, {multAxis, ptAxis}); - histos.add("weight/secMat_mult", "", kTProfile2D, {multAxis, ptAxis}); } } @@ -153,21 +147,48 @@ std::tuple ParticleCompositionCorrection::getWeights(aod::M { static const std::tuple noWeights = {1.f, 1.f, 1.f}; - if (skipAll || std::abs(particle.eta()) > etaCut || particle.pt() < ptMinCut || particle.pt() > ptMaxCut) { + static constexpr int kXi0PDG = 3322; + + /* if (skipAll || std::abs(particle.eta()) > etaCut || particle.pt() < ptMinCut || particle.pt() > ptMaxCut) { + return noWeights; + } */ + + if (skipAll) { return noWeights; } - if (particle.producedByGenerator()) { - if (skipNonPhysicalPrim && !particle.isPhysicalPrimary()) { - return noWeights; - } + if (!particle.has_mothers()) { + return noWeights; + } + + auto particleMother = particle.mothersIds()[0] - particles.offset(); + + // if (particle.producedByGenerator()) { + + if (particleMother == 0. || particleMother == 1. || + (particle.isPhysicalPrimary() && skipNonPhysicalPrim)) { + auto absPDGCode = std::abs(particle.pdgCode()); - // translate abs PDG code to PID variable of neural networks (0: pion, 1: kaon, 2: proton, 3: sigma) if (absPDGCode == PDG_t::kPiPlus || absPDGCode == PDG_t::kPi0) { + + return noWeights; + } + + if (skipNonPhysicalPrim && !particle.isPhysicalPrimary()) { return noWeights; } + if (!particle.has_daughters()) { + + if (particle.eta() > etaCut || particle.pt() > ptMaxCut || particle.pt() < ptMinCut) { + + return noWeights; + } + } + + // translate abs PDG code to PID variable of neural networks (0: pion, 1: kaon, 2: proton, 3: sigma) + static const std::map mapPID = { {PDG_t::kPiPlus, 0.f}, {PDG_t::kPi0, 0.f}, @@ -181,17 +202,17 @@ std::tuple ParticleCompositionCorrection::getWeights(aod::M {PDG_t::kLambda0, 3.f}, {PDG_t::kSigma0, 3.f}, {PDG_t::kXiMinus, 3.f}, - // TODO: potentially extend by xi0/eta/omega/rho/phi/Delta... + {kXi0PDG, 3.f}, // Xi0 + {PDG_t::kOmegaMinus, 3.f}, + // TODO: potentially extend by eta/omega/rho/phi/Delta... // pdg codes defined in AliceO2/Common/Constants/include/CommonConstants/PhysicsConstants.h // e.g. o2::constants::physics::Pdg::kEta }; if (auto iterMapPID = mapPID.find(absPDGCode); iterMapPID != mapPID.end()) { - // LOGP(info, "scaling a {} with status code {} from process {}", particle.pdgCode(), particle.getGenStatusCode(), particle.getProcess()); - float pt = particle.pt(); // calculate particle fractions and corresponding weight for given reference particle - std::vector> input = {{dNdEta}, {pt}, {iterMapPID->second}}; + std::vector> input = {{dNdEta}, {particle.pt()}, {iterMapPID->second}}; float fracData = particleFractionsData.evalModel(input)[0]; float fracMC = particleFractionsMC.evalModel(input)[0]; float weight = (fracMC) ? fracData / fracMC : 1.f; @@ -199,61 +220,88 @@ std::tuple ParticleCompositionCorrection::getWeights(aod::M if (!skipSec && particle.has_daughters()) { storedWeights[particle.index()] = weights; } - if (enableQAHistos && particle.isPhysicalPrimary() && std::abs(particle.eta()) < 0.8) { // o2-linter: disable=magic-number (usual range of charged-partilce measurements) + if (enableQAHistos && particle.isPhysicalPrimary() && std::abs(particle.eta()) < etaCut) { if (iterMapPID->first == PDG_t::kKPlus) { - histos.fill(HIST("frac/data/kaon"), pt, fracData); - histos.fill(HIST("frac/mc/kaon"), pt, fracMC); - histos.fill(HIST("frac/data/kaon_mult"), dNdEta, pt, fracData); - histos.fill(HIST("frac/mc/kaon_mult"), dNdEta, pt, fracMC); - histos.fill(HIST("weight/kaon"), pt, weight); - histos.fill(HIST("weight/kaon_mult"), dNdEta, pt, weight); + histos.fill(HIST("frac/data/kaon"), particle.pt(), fracData); + histos.fill(HIST("frac/mc/kaon"), particle.pt(), fracMC); + histos.fill(HIST("weight/kaon"), particle.pt(), weight); } if (iterMapPID->first == PDG_t::kProton) { - histos.fill(HIST("frac/data/proton"), pt, fracData); - histos.fill(HIST("frac/mc/proton"), pt, fracMC); - histos.fill(HIST("frac/data/proton_mult"), dNdEta, pt, fracData); - histos.fill(HIST("frac/mc/proton_mult"), dNdEta, pt, fracMC); - histos.fill(HIST("weight/proton"), pt, weight); - histos.fill(HIST("weight/proton_mult"), dNdEta, pt, weight); + histos.fill(HIST("frac/data/proton"), particle.pt(), fracData); + histos.fill(HIST("frac/mc/proton"), particle.pt(), fracMC); + histos.fill(HIST("weight/proton"), particle.pt(), weight); } if (iterMapPID->first == PDG_t::kSigmaPlus || iterMapPID->first == PDG_t::kSigmaMinus) { - histos.fill(HIST("frac/data/sigma"), pt, fracData); - histos.fill(HIST("frac/mc/sigma"), pt, fracMC); - histos.fill(HIST("frac/data/sigma_mult"), dNdEta, pt, fracData); - histos.fill(HIST("frac/mc/sigma_mult"), dNdEta, pt, fracMC); - histos.fill(HIST("weight/sigma"), pt, weight); - histos.fill(HIST("weight/sigma_mult"), dNdEta, pt, weight); + histos.fill(HIST("frac/data/sigma"), particle.pt(), fracData); + histos.fill(HIST("frac/mc/sigma"), particle.pt(), fracMC); + histos.fill(HIST("weight/sigma"), particle.pt(), weight); } + + histos.fill(HIST("weight/prim"), particle.pt(), weight); } + return weights; } + } else if (!skipSec) { + auto refParticleID = particle.index(); - // LOGP(error, "Particle [{}] {} from process {}", refParticleID, particle.pdgCode(), particle.getProcess()); - while (!particles.iteratorAt(refParticleID).producedByGenerator() && particles.iteratorAt(refParticleID).has_mothers()) { + auto pdgParticle = pdg->GetParticle(particle.pdgCode()); + + if (!pdgParticle) { + // LOGP(info, "Unknown PDG code: {}", particle.pdgCode()); + // return noWeights; + } else { + + if (skipNeutralSec && pdgParticle->Charge() == 0.) { + return noWeights; + } + } + + if (particle.eta() > etaCut || particle.pt() > ptMaxCut || particle.pt() < ptMinCut) { + return noWeights; + } + + while (particles.iteratorAt(refParticleID).has_mothers()) { + auto motherID = particles.iteratorAt(refParticleID).mothersIds()[0] - particles.offset(); - // LOGP(error, "-> mom [{}] {} from process {}", motherID, particles.iteratorAt(motherID).pdgCode(), particles.iteratorAt(motherID).getProcess()); + + auto mom = particles.iteratorAt(motherID); + + if (!mom.has_mothers()) { + break; + } + refParticleID = motherID; + + if (mom.isPhysicalPrimary() && skipNonPhysicalPrim) { + break; + } + } + + if (refParticleID == 0 || refParticleID == 1) { + LOGF(error, "GLOBAL MOTHER IS INCIDING PROTON"); } if (storedWeights.find(refParticleID) == storedWeights.end()) { - // LOGP(error, " no ref particle stored for particle {} from process {}!!", particle.pdgCode(), particle.getProcess()); + // LOGP(info, "MCB: no ref stored for mom {}, PDG: {}!!", particles.iteratorAt(refParticleID).index(),particles.iteratorAt(refParticleID).pdgCode()); return noWeights; } if (enableQAHistos) { + float weight = get<0>(storedWeights.at(refParticleID)); - auto pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (pdgParticle && pdgParticle->Charge() != 0.) { - if (particle.getProcess() == TMCProcess::kPDecay) { + if (particle.getProcess() == TMCProcess::kPDecay || particle.getProcess() == TMCProcess::kPPrimary) { histos.fill(HIST("weight/secDec"), particle.pt(), weight); - histos.fill(HIST("weight/secDec_mult"), dNdEta, particle.pt(), weight); - } else if (particle.getProcess() == TMCProcess::kPHInhelastic || particle.getProcess() == TMCProcess::kPHadronic || particle.getProcess() == TMCProcess::kPHElastic) { + // } else if (particle.getProcess() == TMCProcess::kPHInhelastic || particle.getProcess() == TMCProcess::kPHadronic || particle.getProcess() == TMCProcess::kPHElastic) { + } else { histos.fill(HIST("weight/secMat"), particle.pt(), weight); - histos.fill(HIST("weight/secMat_mult"), dNdEta, particle.pt(), weight); } } } + return storedWeights.at(refParticleID); } return noWeights; @@ -266,6 +314,7 @@ void ParticleCompositionCorrection::process(aod::McCollisions::iterator const&, float dNdEtafid = 0.f; for (const auto& particle : particles) { + if (!particle.isPhysicalPrimary()) { continue; } @@ -288,6 +337,8 @@ void ParticleCompositionCorrection::process(aod::McCollisions::iterator const&, histos.fill(HIST("multDist_fidVsINELg0"), dNdEta, dNdEtafid); } + // LOGP(info, "MCB: ------- New event with dNdeta = {} in |eta| < 0.8 -------", dNdEtafid); + std::map> storedWeights; for (const auto& particle : particles) { auto [weight, weightSysUp, weightSysDown] = getWeights(particles, particle, storedWeights, dNdEta); diff --git a/PWGLF/Tasks/Nuspex/pccQa.cxx b/PWGLF/Tasks/Nuspex/pccQa.cxx index fe04b0ffc4e..f2843e5127d 100644 --- a/PWGLF/Tasks/Nuspex/pccQa.cxx +++ b/PWGLF/Tasks/Nuspex/pccQa.cxx @@ -30,18 +30,76 @@ #include #include +#include #include +#include #include using namespace o2; using namespace o2::framework; using aod::track::TrackSelectionFlags; +int getMotherCategory(int pdg) +{ + switch (std::abs(pdg)) { + case 310: + return 1; // K0S + case 130: + return 2; // K0L + case 3122: + return 3; // Lambda + + case 3222: // Sigma+ + case 3112: + return 4; // Sigma- + + case 3312: + return 5; // Xi- + + case 3334: + return 6; // Omega + + case 211: + return 7; // Pion + case 321: + return 8; // Kaon + case 2112: + return 9; // Neutron + case 13: + return 10; // Muon + + default: + return 11; // Other + } +} + +static std::string getMotherName(int category) +{ + static const std::vector names = { + "K0S", + "K0L", + "Lambda", + "SigmaCharged", + "XiCharged", + "Omega", + "Pi", + "K", + "n", + "Muon", + "Other"}; + + return names.at(category - 1); +} + struct PccQa { HistogramRegistry histos; Service pdg; + std::vector> motherDCA; + std::vector> motherDCAWeighted; + std::vector> motherWeights; + static constexpr float MaxVtxZ = 10.f; Configurable maxMult{"maxMult", 100, "max multiplicity"}; @@ -73,44 +131,76 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) void PccQa::init(InitContext const&) { - histos.add("eventCounter", "", kTH1D, {{1, 0.5, 1.5}}); - const AxisSpec dcaAxis = {1000, -1., 1., "#it{DCA}_{xy}", "dca"}; + // const AxisSpec dcaAxis = {1000, -1., 1., "#it{DCA}_{xy}", "dca"}; + const AxisSpec dcaAxis = {99, -1., 1., "#it{DCA}_{xy}", "dca"}; std::vector ptBinEdges = {0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0}; const AxisSpec ptAxis{ptBinEdges, "#it{p}_{T} (GeV/#it{c})", "pt"}; - const int nBinsMult = 30; - const AxisSpec nchAxis = {nBinsMult, 0., static_cast(maxMult), "#it{N}_{ch}", "nch"}; + const int nBinsMult = 25; + const AxisSpec nchAxis = {nBinsMult, 1., static_cast(maxMult), "#it{N}_{ch}", "nch"}; + + histos.add("multDist", "", kTH1D, {nchAxis}); - histos.add("DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis}); histos.add("DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); if (doprocessMC) { - histos.add("DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis}); histos.add("DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); - histos.add("prim/DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis}); histos.add("prim/DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); - histos.add("prim/DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); - histos.add("prim/DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis}); - histos.add("sec/DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis}); histos.add("sec/DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); - - histos.add("sec/DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis}); histos.add("sec/DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); - histos.add("sec/dec/DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis}); histos.add("sec/dec/DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); - - histos.add("sec/dec/DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis}); histos.add("sec/dec/DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); - histos.add("sec/mat/DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis}); histos.add("sec/mat/DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); - - histos.add("sec/mat/DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis}); histos.add("sec/mat/DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis}); + + histos.add("prim/weights", "", kTProfile, {ptAxis}); + histos.add("sec/weights", "", kTProfile, {ptAxis}); + histos.add("sec/dec/weights", "", kTProfile, {ptAxis}); + histos.add("sec/mat/weights", "", kTProfile, {ptAxis}); + + std::vector labels = { + "K^{0}_{S}", "K^{0}_{L}", "#Lambda", + "#Sigma^{#pm}", + "#Xi^{#pm}", + "#Omega^{#pm}", + "#pi^{#pm}", "K^{#pm}", "n", + "#mu^{#pm}", + "Other"}; + + const int nMotherCategories = static_cast(labels.size()); + + histos.add( + "sec/dec/motherPDG", + "Mother particles of decay secondaries", + kTH1D, + {{nMotherCategories, 0.5, nMotherCategories + 0.5}}); + + auto hMother = histos.get(HIST("sec/dec/motherPDG")); + + for (int i = 0; i < nMotherCategories; ++i) { + hMother->GetXaxis()->SetBinLabel(i + 1, labels[i].c_str()); + } + + for (int category = 1; category <= nMotherCategories; ++category) { + const auto mother = getMotherName(category); + + motherDCA.push_back(histos.add( + fmt::format("sec/dec/{}/DCAxyVsPt", mother), + "", kTH2D, {ptAxis, dcaAxis})); + + motherDCAWeighted.push_back(histos.add( + fmt::format("sec/dec/{}/DCAxyVsPt_weighted", mother), + "", kTH2D, {ptAxis, dcaAxis})); + + motherWeights.push_back(histos.add( + fmt::format("sec/dec/{}/weights", mother), + "", kTProfile, {ptAxis})); + } } } @@ -135,44 +225,109 @@ void PccQa::processMeas(const C& collision, const T& tracks) if ((std::abs(collision.posZ()) > MaxVtxZ) || !collision.sel8()) { return; } - histos.fill(HIST("eventCounter"), 1); + + auto MultCounterMeas = 0; for (const auto& track : tracks) { + if (!TrackSelectionFlags::checkFlag(track.trackCutFlag(), TrackSelectionFlags::kGlobalTrackWoDCA)) { continue; } - histos.fill(HIST("DCAxyVsPt"), track.pt(), track.dcaXY()); - histos.fill(HIST("DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY()); + + MultCounterMeas++; + } + + /* LOGP(info, "MCB: --------------------------"); + LOGP(info, "MCB: ++ New event with Nch = {} ",MultCounterMeas); + LOGP(info, "MCB: --------------------------"); */ + + histos.fill(HIST("multDist"), MultCounterMeas); + + for (const auto& track : tracks) { + if (!TrackSelectionFlags::checkFlag(track.trackCutFlag(), TrackSelectionFlags::kGlobalTrackWoDCA)) { + continue; + } + + histos.fill(HIST("DCAxyVsPtVsMult"), MultCounterMeas, track.pt(), track.dcaXY()); if constexpr (IS_MC) { if (!track.has_mcParticle()) { continue; } + const auto& particle = track.template mcParticle_as(); - histos.fill(HIST("DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight()); - histos.fill(HIST("DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight()); + histos.fill(HIST("DCAxyVsPtVsMult_weighted"), MultCounterMeas, track.pt(), track.dcaXY(), particle.pccWeight()); if (particle.isPhysicalPrimary()) { - histos.fill(HIST("prim/DCAxyVsPt"), track.pt(), track.dcaXY()); - histos.fill(HIST("prim/DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight()); - histos.fill(HIST("prim/DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY()); - histos.fill(HIST("prim/DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight()); + + histos.fill(HIST("prim/weights"), particle.pt(), particle.pccWeight()); + histos.fill(HIST("prim/DCAxyVsPtVsMult"), MultCounterMeas, track.pt(), track.dcaXY()); + histos.fill(HIST("prim/DCAxyVsPtVsMult_weighted"), MultCounterMeas, track.pt(), track.dcaXY(), particle.pccWeight()); + } else { - histos.fill(HIST("sec/DCAxyVsPt"), track.pt(), track.dcaXY()); - histos.fill(HIST("sec/DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight()); - histos.fill(HIST("sec/DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY()); - histos.fill(HIST("sec/DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight()); - if (particle.getProcess() == TMCProcess::kPDecay) { - histos.fill(HIST("sec/dec/DCAxyVsPt"), track.pt(), track.dcaXY()); - histos.fill(HIST("sec/dec/DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight()); - histos.fill(HIST("sec/dec/DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY()); - histos.fill(HIST("sec/dec/DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight()); - } else if (particle.getProcess() == TMCProcess::kPHInhelastic || particle.getProcess() == TMCProcess::kPHadronic || particle.getProcess() == TMCProcess::kPHElastic) { - histos.fill(HIST("sec/mat/DCAxyVsPt"), track.pt(), track.dcaXY()); - histos.fill(HIST("sec/mat/DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight()); - histos.fill(HIST("sec/mat/DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY()); - histos.fill(HIST("sec/mat/DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight()); + + if (!particle.has_mothers()) { + continue; + } + + histos.fill(HIST("sec/weights"), particle.pt(), particle.pccWeight()); + + histos.fill(HIST("sec/DCAxyVsPtVsMult"), MultCounterMeas, track.pt(), track.dcaXY()); + histos.fill(HIST("sec/DCAxyVsPtVsMult_weighted"), MultCounterMeas, track.pt(), track.dcaXY(), particle.pccWeight()); + + if (particle.getProcess() == TMCProcess::kPDecay || particle.getProcess() == TMCProcess::kPPrimary) { + + histos.fill(HIST("sec/dec/DCAxyVsPtVsMult"), + MultCounterMeas, track.pt(), track.dcaXY()); + + histos.fill(HIST("sec/dec/DCAxyVsPtVsMult_weighted"), + MultCounterMeas, track.pt(), track.dcaXY(), + particle.pccWeight()); + + histos.fill(HIST("sec/dec/weights"), + particle.pt(), particle.pccWeight()); + + /* LOGP(info, + " -> Particle={}, PDG={}, process={}, " + "nMoms={}, 1. mom ID={}", + particle.globalIndex(), + particle.pdgCode(), + particle.getProcess(), + particle.mothersIds().size(), + particle.mothersIds()[0]); */ + + // Identify the immediate mother + // auto mother = particle.template mothers_first_as(); + auto mother = particle.template mothers_first_as(); + + const int motherCategory = getMotherCategory(mother.pdgCode()); + const std::string motherName = getMotherName(motherCategory); + + histos.fill( + HIST("sec/dec/motherPDG"), + motherCategory); + + const int i = motherCategory - 1; + + /* LOGP(info, + " --> Mother PDG={}", + mother.pdgCode() + ); */ + + motherDCA[i]->Fill(track.pt(), track.dcaXY()); + + motherDCAWeighted[i]->Fill( + track.pt(), track.dcaXY(), particle.pccWeight()); + + motherWeights[i]->Fill( + particle.pt(), particle.pccWeight()); + + } else { + + histos.fill(HIST("sec/mat/weights"), particle.pt(), particle.pccWeight()); + histos.fill(HIST("sec/mat/DCAxyVsPtVsMult"), MultCounterMeas, track.pt(), track.dcaXY()); + histos.fill(HIST("sec/mat/DCAxyVsPtVsMult_weighted"), MultCounterMeas, track.pt(), track.dcaXY(), particle.pccWeight()); } } }