diff --git a/ALICE3/Tasks/alice3-qa-singleparticle.cxx b/ALICE3/Tasks/alice3-qa-singleparticle.cxx index 25d21fcc7cd..851f67adbe7 100644 --- a/ALICE3/Tasks/alice3-qa-singleparticle.cxx +++ b/ALICE3/Tasks/alice3-qa-singleparticle.cxx @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// -/// \file qa-singleparticle.cxx +/// \file alice3-qa-singleparticle.cxx /// \author Nicolò Jacazio nicolo.jacazio@cern.ch /// \brief Task to monitor the single particle QA, at the particle and track level, showing the tracked and the origin of particles /// @@ -28,7 +28,6 @@ #include #include -#include #include #include #include @@ -38,7 +37,9 @@ #include #include +#include #include +#include #include using namespace o2; @@ -46,26 +47,50 @@ using namespace o2::framework; using namespace o2::framework::expressions; struct Alice3SingleParticle { - Service pdg; - Configurable PDG{"PDG", 2212, "PDG code of the particle of interest"}; - Configurable IsStable{"IsStable", 0, "Flag to check stable particles"}; + Service pdg{}; HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - Configurable ptBins{"pt-bins", 500, "Number of pT bins"}; - Configurable ptMin{"pt-min", 0.f, "Lower limit in pT"}; - Configurable ptMax{"pt-max", 10.f, "Upper limit in pT"}; - Configurable etaBins{"eta-bins", 500, "Number of eta bins"}; - Configurable etaMin{"eta-min", -3.f, "Lower limit in eta"}; - Configurable etaMax{"eta-max", 3.f, "Upper limit in eta"}; - Configurable yMin{"y-min", -3.f, "Lower limit in y"}; - Configurable yMax{"y-max", 3.f, "Upper limit in y"}; - Configurable prodBins{"prod-bins", 100, "Number of production vertex bins"}; - Configurable prodMin{"prod-min", -1.f, "Lower limit in production vertex"}; - Configurable prodMax{"prod-max", 1.f, "Upper limit in production vertex"}; - Configurable prodBinsZ{"prod-bins-z", 100, "Number of production vertex bins"}; - Configurable prodMinZ{"prod-min-z", -10.f, "Lower limit in production vertex along Z"}; - Configurable prodMaxZ{"prod-max-z", 10.f, "Upper limit in production vertex along Z"}; - Configurable charge{"charge", 1.f, "Particle charge to scale the reconstructed momentum"}; - Configurable doPrint{"doPrint", false, "Flag to print debug messages"}; + Configurable pdgCode{"pdgCode", 2212, "PDG code of the particle of interest"}; + + struct : ConfigurableGroup { + std::string prefix = "settings"; + Configurable isStable{"isStable", 0, "Flag to check stable particles"}; + Configurable charge{"charge", 1.f, "Particle charge to scale the reconstructed momentum"}; + Configurable nchEtaMax{"nchEtaMax", 0.5f, "|eta| acceptance of the Nch estimator (use the Ditto tune mActivityEtaMax)"}; + Configurable doPrint{"doPrint", false, "Flag to print debug messages"}; + } settings; + + // The pT, eta, y and production-vertex-z axes also define the selection of the particles of interest + struct : ConfigurableGroup { + std::string prefix = "axes"; + ConfigurableAxis axisMomentum{"axisMomentum", {500, 0.f, 10.f}, "Binning for p and pT"}; + ConfigurableAxis axisMomentumComponent{"axisMomentumComponent", {500, -10.f, 10.f}, "Binning for px, py and pz"}; + ConfigurableAxis axisEnergy{"axisEnergy", {500, 0.f, 1000.f}, "Binning for the energy"}; + ConfigurableAxis axisEta{"axisEta", {500, -3.f, 3.f}, "Binning for eta"}; + ConfigurableAxis axisY{"axisY", {200, -3.f, 3.f}, "Binning for rapidity"}; + ConfigurableAxis axisVtxXY{"axisVtxXY", {100, -1.f, 1.f}, "Binning for the event vertex in X and Y"}; + ConfigurableAxis axisVtxZ{"axisVtxZ", {100, -20.f, 20.f}, "Binning for the event vertex in Z"}; + ConfigurableAxis axisProdXY{"axisProdXY", {100, -1.f, 1.f}, "Binning for the production vertex in X and Y"}; + ConfigurableAxis axisProdZ{"axisProdZ", {100, -10.f, 10.f}, "Binning for the production vertex in Z"}; + ConfigurableAxis axisProdRadius{"axisProdRadius", {100, 0.f, 2.f}, "Binning for the production vertex radius"}; + ConfigurableAxis axisNch{"axisNch", {150, 0.f, 150.f}, "Charged-particle multiplicity"}; + ConfigurableAxis axisNSignal{"axisNSignal", {20, -0.5f, 19.5f}, "Number of selected particles of interest per collision"}; + } axes; + + std::vector particlesOfInterest; // global indices of the MC particles selected in genericProcessParticles, per DF + + // Selection range [min, max) of the particles of interest, taken from an axis so that every selected particle is inside its histograms + struct Range { + double min = 0.; + double max = 0.; + bool contains(double value) const { return value >= min && value < max; } + }; + Range ptRange, etaRange, yRange, prodZRange; + + static Range getRange(const ConfigurableAxis& axis) + { + const AxisSpec spec{axis}; + return {spec.binEdges.front(), spec.binEdges.back()}; + } void init(InitContext&) { @@ -95,123 +120,123 @@ struct Alice3SingleParticle { pdg->AddParticle("helium4", "helium4", o2::constants::physics::MassAlpha, kTRUE, 0.0, 6, "Nucleus", 1000020040); pdg->AddAntiParticle("anti-helium4", -1000020040); - const TString tit = Form("%i", PDG.value); - const AxisSpec axisVx{100, -1, 1, "Vtx_{x}"}; - const AxisSpec axisVy{100, -1, 1, "Vtx_{y}"}; - const AxisSpec axisVz{100, -20, 20, "Vtx_{z}"}; + ptRange = getRange(axes.axisMomentum); + etaRange = getRange(axes.axisEta); + yRange = getRange(axes.axisY); + prodZRange = getRange(axes.axisProdZ); + LOG(info) << "Selecting particles of interest in pT [" << ptRange.min << ", " << ptRange.max << "), eta [" << etaRange.min << ", " << etaRange.max + << "), y [" << yRange.min << ", " << yRange.max << "), production z [" << prodZRange.min << ", " << prodZRange.max << ")"; + + const TString titlePdgCode = Form("%i", pdgCode.value); const AxisSpec axisPDGs{100, 0.f, 100.f, "PDG code"}; const AxisSpec axisCharge{21, -10.5f, 10.5f, "Charge (e/3)"}; - const AxisSpec axisP{ptBins, ptMax, ptMax, "#it{p} (GeV/#it{c})"}; - const AxisSpec axisPt{ptBins, ptMin, ptMax, "#it{p}_{T} (GeV/#it{c})"}; - const AxisSpec axisPtMC{ptBins, ptMin, ptMax, "#it{p}_{T}^{MC} (GeV/#it{c})"}; - const AxisSpec axisPx{ptBins, -ptMax, ptMax, "#it{p}_{x} (GeV/#it{c})"}; - const AxisSpec axisPy{ptBins, -ptMax, ptMax, "#it{p}_{y} (GeV/#it{c})"}; - const AxisSpec axisPz{ptBins, -ptMax, ptMax, "#it{p}_{z} (GeV/#it{c})"}; - const AxisSpec axisEta{etaBins, etaMin, etaMax, "#it{#eta}"}; - const AxisSpec axisY{etaBins, yMin, yMax, "#it{y}"}; - const AxisSpec axisE{etaBins, 0, 1000, "E"}; - const AxisSpec axisProdx{prodBins, prodMin, prodMax, "Prod. Vertex X (cm)"}; - const AxisSpec axisPrody{prodBins, prodMin, prodMax, "Prod. Vertex Y (cm)"}; - const AxisSpec axisProdz{prodBinsZ, prodMinZ, prodMaxZ, "Prod. Vertex Z (cm)"}; - const AxisSpec axisProdRadius{prodBins, 0., 2. * prodMax, "Prod. Vertex Radius (cm)"}; - - histos.add("event/VtxX", "Vertex X", kTH1D, {axisVx}); - histos.add("event/VtxY", "Vertex Y", kTH1D, {axisVy}); - histos.add("event/VtxZ", "Vertex Z", kTH1D, {axisVz}); + + histos.add("AnalysedParticlePdgCode", "Analysed particle PDG code", kTH1D, {{1, -0.5, 0.5}}); + histos.get(HIST("AnalysedParticlePdgCode"))->SetBinContent(1, pdgCode); + + const auto* particle = pdg->GetParticle(pdgCode); + histos.add("AnalysedParticleMass", "Analysed particle mass", kTH1D, {{1, -0.5, 0.5}}); + histos.get(HIST("AnalysedParticleMass"))->SetBinContent(1, particle ? particle->Mass() : -1.); + + histos.add("event/VtxX", "Vertex X", kTH1D, {axes.axisVtxXY}); + histos.add("event/VtxY", "Vertex Y", kTH1D, {axes.axisVtxXY}); + histos.add("event/VtxZ", "Vertex Z", kTH1D, {axes.axisVtxZ}); + histos.add("event/Nch", "Events;N_{ch};Events", kTH1D, {axes.axisNch}); + histos.add("particle/NSignalVsNch", "Particles of interest per event " + titlePdgCode + ";N_{ch};N_{signal}", kTH2D, {axes.axisNch, axes.axisNSignal}); histos.add("particle/PDGs", "Particle PDGs", kTH2D, {axisPDGs, axisCharge}); histos.add("particle/PDGsPrimaries", "Particle PDGs of Primaries", kTH2D, {axisPDGs, axisCharge}); histos.add("particle/PDGsSecondaries", "Particle PDGs of Secondaries", kTH2D, {axisPDGs, axisCharge}); - auto h = histos.add("particle/producedProcess", "Particle process of production " + tit, kTH1D, {{kMaxMCProcess, -0.5, -0.5 + kMaxMCProcess, "getProcess"}}); + auto h = histos.add("particle/producedProcess", "Particle process of production " + titlePdgCode, kTH1D, {{kMaxMCProcess, -0.5, -0.5 + kMaxMCProcess, "getProcess"}}); for (int i = 0; i < kMaxMCProcess; i++) { h->GetXaxis()->SetBinLabel(i + 1, TMCProcessName[i]); } - histos.add("particle/Pt", "Particle Pt " + tit, kTH1D, {axisPt}); - histos.add("particle/P", "Particle P " + tit, kTH1D, {axisP}); - histos.add("particle/primariesPt", "Particle Pt (primary) " + tit, kTH1D, {axisPt}); - histos.add("particle/primariesP", "Particle P (primary) " + tit, kTH1D, {axisP}); - histos.add("particle/secondariesPt", "Particle Pt (secondary) " + tit, kTH1D, {axisPt}); - histos.add("particle/secondariesP", "Particle P (secondary) " + tit, kTH1D, {axisP}); - histos.add("particle/prodVx", "Particle Prod. Vertex X " + tit, kTH1D, {axisProdx}); - histos.add("particle/prodVy", "Particle Prod. Vertex Y " + tit, kTH1D, {axisPrody}); - histos.add("particle/prodVz", "Particle Prod. Vertex Z " + tit, kTH1D, {axisProdz}); - histos.add("particle/prodRadius", "Particle Prod. Vertex Radius " + tit, kTH1D, {axisProdRadius}); - histos.add("particle/prodVxVy", "Particle Prod. Vertex X vs Vertex Y " + tit, kTH2D, {axisProdx, axisPrody}); - histos.add("particle/prodVxVyPrm", "Particle Prod. Vertex X vs Vertex Y Primary " + tit, kTH2D, {axisProdx, axisPrody}); - histos.add("particle/prodVxVyStr", "Particle Prod. Vertex X vs Vertex Y Decay " + tit, kTH2D, {axisProdx, axisPrody}); - histos.add("particle/prodVxVyMat", "Particle Prod. Vertex X vs Vertex Y Material " + tit, kTH2D, {axisProdx, axisPrody}); - histos.add("particle/prodVxVsPt", "Particle Prod. Vertex X " + tit, kTH2D, {axisPt, axisProdx}); - histos.add("particle/prodVyVsPt", "Particle Prod. Vertex Y " + tit, kTH2D, {axisPt, axisPrody}); - histos.add("particle/prodVzVsPt", "Particle Prod. Vertex Z " + tit, kTH2D, {axisPt, axisProdz}); - histos.add("particle/prodRadiusVsPt", "Particle Prod. Vertex Radius " + tit, kTH2D, {axisPt, axisProdRadius}); - histos.add("particle/prodRadiusVsEta", "Particle Prod. Vertex Radius " + tit, kTH2D, {axisEta, axisProdRadius}); - histos.add("particle/prodRadius3DVsPt", "Particle Prod. Vertex Radius XYZ " + tit, kTH2D, {axisPt, axisProdRadius}); - histos.add("particle/Eta", "Particle Eta " + tit, kTH1D, {axisEta}); - histos.add("particle/primariesEta", "Particle Eta (primary) " + tit, kTH1D, {axisEta}); - histos.add("particle/secondariesEta", "Particle Eta (secondary) " + tit, kTH1D, {axisEta}); - histos.add("particle/Y", "Particle Y " + tit, kTH1D, {axisY}); - histos.add("particle/primariesY", "Particle Y (primary)" + tit, kTH1D, {axisY}); - histos.add("particle/secondariesY", "Particle Y (secondary)" + tit, kTH1D, {axisY}); - histos.add("particle/EvsPz", "Particle E vs Pz " + tit, kTH2D, {axisE, axisPz}); - histos.add("particle/YvzPz", "Particle Y vs Pz " + tit, kTH2D, {axisY, axisPz}); - histos.add("particle/EtavzPz", "Particle Eta vs Pz " + tit, kTH2D, {axisEta, axisPz}); - histos.add("particle/PtvzPz", "Particle Pt vs Pz " + tit, kTH2D, {axisPt, axisPz}); - histos.add("particle/PvzPz", "Particle P vs Pz " + tit, kTH2D, {axisP, axisPz}); - histos.add("particle/Px", "Particle Px " + tit, kTH1D, {axisPx}); - histos.add("particle/Py", "Particle Py " + tit, kTH1D, {axisPy}); - histos.add("particle/Pz", "Particle Pz " + tit, kTH1D, {axisPz}); - histos.add("particle/decayDistanceX", "Particle decay distance " + tit, kTH1D, {axisProdx}); - histos.add("particle/decayDistanceY", "Particle decay distance " + tit, kTH1D, {axisPrody}); - histos.add("particle/decayDistanceZ", "Particle decay distance " + tit, kTH1D, {axisProdz}); - histos.add("particle/decayDistanceVsPt", "Particle decay distance " + tit, kTH2D, {axisPt, axisProdRadius}); - histos.add("particle/decayCtau", "Particle decay distance " + tit, kTH1D, {{1000, 0, 10000, "c#tau (#mum)"}}); - - histos.add("particle/daughters/Number", "Number of Daughters " + tit, kTH1D, {{20, -0.5, 19.5}}); - histos.add("particle/daughters/PDGs", "Daughters PDGs " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("particle/daughters/PDGsPrimaries", "Daughters PDGs Primaries of " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("particle/daughters/PDGsSecondaries", "Daughters PDGs Secondaries of " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("particle/daughters/prodVx", "Daughters Prod. Vertex X " + tit, kTH1D, {axisProdx}); - histos.add("particle/daughters/prodVy", "Daughters Prod. Vertex Y " + tit, kTH1D, {axisPrody}); - histos.add("particle/daughters/prodVz", "Daughters Prod. Vertex Z " + tit, kTH1D, {axisProdz}); - histos.add("particle/daughters/prodRadiusVsPt", "Daughters Prod. Vertex Radius " + tit, kTH2D, {axisPt, axisProdRadius}); - histos.add("particle/daughters/prodRadius3DVsPt", "Daughters Prod. Vertex Radius XYZ " + tit, kTH2D, {axisPt, axisProdRadius}); - - histos.add("particle/mothers/Number", "Number of Mothers " + tit, kTH1D, {{20, -0.5, 19.5}}); - histos.add("particle/mothers/PDGs", "Mothers PDGs " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("particle/mothers/PDGsPrimaries", "Mothers PDGs Primaries of " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("particle/mothers/PDGsSecondaries", "Mothers PDGs Secondaries of " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("particle/mothers/prodVx", "Mothers Prod. Vertex X " + tit, kTH1D, {axisProdx}); - histos.add("particle/mothers/prodVy", "Mothers Prod. Vertex Y " + tit, kTH1D, {axisPrody}); - histos.add("particle/mothers/prodVz", "Mothers Prod. Vertex Z " + tit, kTH1D, {axisProdz}); - histos.add("particle/mothers/prodRadiusVsPt", "Mothers Prod. Vertex Radius " + tit, kTH2D, {axisPt, axisProdRadius}); - histos.add("particle/mothers/prodRadius3DVsPt", "Mothers Prod. Vertex Radius XYZ " + tit, kTH2D, {axisPt, axisProdRadius}); + histos.add("particle/Pt", "Particle Pt " + titlePdgCode, kTH1D, {axes.axisMomentum}); + histos.add("particle/P", "Particle P " + titlePdgCode, kTH1D, {axes.axisMomentum}); + histos.add("particle/primariesPt", "Particle Pt (primary) " + titlePdgCode, kTH1D, {axes.axisMomentum}); + histos.add("particle/primariesP", "Particle P (primary) " + titlePdgCode, kTH1D, {axes.axisMomentum}); + histos.add("particle/secondariesPt", "Particle Pt (secondary) " + titlePdgCode, kTH1D, {axes.axisMomentum}); + histos.add("particle/secondariesP", "Particle P (secondary) " + titlePdgCode, kTH1D, {axes.axisMomentum}); + histos.add("particle/prodVx", "Particle Prod. Vertex X " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/prodVy", "Particle Prod. Vertex Y " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/prodVz", "Particle Prod. Vertex Z " + titlePdgCode, kTH1D, {axes.axisProdZ}); + histos.add("particle/prodRadius", "Particle Prod. Vertex Radius " + titlePdgCode, kTH1D, {axes.axisProdRadius}); + histos.add("particle/prodVxVy", "Particle Prod. Vertex X vs Vertex Y " + titlePdgCode, kTH2D, {axes.axisProdXY, axes.axisProdXY}); + histos.add("particle/prodVxVyPrm", "Particle Prod. Vertex X vs Vertex Y Primary " + titlePdgCode, kTH2D, {axes.axisProdXY, axes.axisProdXY}); + histos.add("particle/prodVxVyStr", "Particle Prod. Vertex X vs Vertex Y Decay " + titlePdgCode, kTH2D, {axes.axisProdXY, axes.axisProdXY}); + histos.add("particle/prodVxVyMat", "Particle Prod. Vertex X vs Vertex Y Material " + titlePdgCode, kTH2D, {axes.axisProdXY, axes.axisProdXY}); + histos.add("particle/prodVxVsPt", "Particle Prod. Vertex X " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdXY}); + histos.add("particle/prodVyVsPt", "Particle Prod. Vertex Y " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdXY}); + histos.add("particle/prodVzVsPt", "Particle Prod. Vertex Z " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdZ}); + histos.add("particle/prodRadiusVsPt", "Particle Prod. Vertex Radius " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdRadius}); + histos.add("particle/prodRadiusVsEta", "Particle Prod. Vertex Radius " + titlePdgCode, kTH2D, {axes.axisEta, axes.axisProdRadius}); + histos.add("particle/prodRadius3DVsPt", "Particle Prod. Vertex Radius XYZ " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdRadius}); + histos.add("particle/Eta", "Particle Eta " + titlePdgCode, kTH1D, {axes.axisEta}); + histos.add("particle/primariesEta", "Particle Eta (primary) " + titlePdgCode, kTH1D, {axes.axisEta}); + histos.add("particle/secondariesEta", "Particle Eta (secondary) " + titlePdgCode, kTH1D, {axes.axisEta}); + histos.add("particle/Y", "Particle Y " + titlePdgCode, kTH1D, {axes.axisY}); + histos.add("particle/primariesY", "Particle Y (primary)" + titlePdgCode, kTH1D, {axes.axisY}); + histos.add("particle/secondariesY", "Particle Y (secondary)" + titlePdgCode, kTH1D, {axes.axisY}); + histos.add("particle/EvsPz", "Particle E vs Pz " + titlePdgCode, kTH2D, {axes.axisEnergy, axes.axisMomentumComponent}); + histos.add("particle/YvzPz", "Particle Y vs Pz " + titlePdgCode, kTH2D, {axes.axisY, axes.axisMomentumComponent}); + histos.add("particle/EtavzPz", "Particle Eta vs Pz " + titlePdgCode, kTH2D, {axes.axisEta, axes.axisMomentumComponent}); + histos.add("particle/PtvzPz", "Particle Pt vs Pz " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisMomentumComponent}); + histos.add("particle/PvzPz", "Particle P vs Pz " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisMomentumComponent}); + histos.add("particle/Px", "Particle Px " + titlePdgCode, kTH1D, {axes.axisMomentumComponent}); + histos.add("particle/Py", "Particle Py " + titlePdgCode, kTH1D, {axes.axisMomentumComponent}); + histos.add("particle/Pz", "Particle Pz " + titlePdgCode, kTH1D, {axes.axisMomentumComponent}); + histos.add("particle/decayDistanceX", "Particle decay distance " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/decayDistanceY", "Particle decay distance " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/decayDistanceZ", "Particle decay distance " + titlePdgCode, kTH1D, {axes.axisProdZ}); + histos.add("particle/decayDistanceVsPt", "Particle decay distance " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdRadius}); + histos.add("particle/decayCtau", "Particle decay distance " + titlePdgCode, kTH1D, {{1000, 0, 10000, "c#tau (#mum)"}}); + + histos.add("particle/daughters/Number", "Number of Daughters " + titlePdgCode, kTH1D, {{20, -0.5, 19.5}}); + histos.add("particle/daughters/PDGs", "Daughters PDGs " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("particle/daughters/PDGsPrimaries", "Daughters PDGs Primaries of " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("particle/daughters/PDGsSecondaries", "Daughters PDGs Secondaries of " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("particle/daughters/prodVx", "Daughters Prod. Vertex X " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/daughters/prodVy", "Daughters Prod. Vertex Y " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/daughters/prodVz", "Daughters Prod. Vertex Z " + titlePdgCode, kTH1D, {axes.axisProdZ}); + histos.add("particle/daughters/prodRadiusVsPt", "Daughters Prod. Vertex Radius " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdRadius}); + histos.add("particle/daughters/prodRadius3DVsPt", "Daughters Prod. Vertex Radius XYZ " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdRadius}); + + histos.add("particle/mothers/Number", "Number of Mothers " + titlePdgCode, kTH1D, {{20, -0.5, 19.5}}); + histos.add("particle/mothers/PDGs", "Mothers PDGs " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("particle/mothers/PDGsPrimaries", "Mothers PDGs Primaries of " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("particle/mothers/PDGsSecondaries", "Mothers PDGs Secondaries of " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("particle/mothers/prodVx", "Mothers Prod. Vertex X " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/mothers/prodVy", "Mothers Prod. Vertex Y " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/mothers/prodVz", "Mothers Prod. Vertex Z " + titlePdgCode, kTH1D, {axes.axisProdZ}); + histos.add("particle/mothers/prodRadiusVsPt", "Mothers Prod. Vertex Radius " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdRadius}); + histos.add("particle/mothers/prodRadius3DVsPt", "Mothers Prod. Vertex Radius XYZ " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdRadius}); // Go up one generation - histos.add("particle/mothers/mothers/Number", "Number of Mothers mothers " + tit, kTH1D, {{20, -0.5, 19.5}}); - histos.add("particle/mothers/mothers/PDGs", "Mothers mothers PDGs " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("particle/mothers/mothers/PDGsPrimaries", "Mothers mothers PDGs Primaries of " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("particle/mothers/mothers/PDGsSecondaries", "Mothers mothers PDGs Secondaries of " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("particle/mothers/mothers/prodVx", "Mothers mothers Prod. Vertex X " + tit, kTH1D, {axisProdx}); - histos.add("particle/mothers/mothers/prodVy", "Mothers mothers Prod. Vertex Y " + tit, kTH1D, {axisPrody}); - histos.add("particle/mothers/mothers/prodVz", "Mothers mothers Prod. Vertex Z " + tit, kTH1D, {axisProdz}); - histos.add("particle/mothers/mothers/prodRadiusVsPt", "Mothers mothers Prod. Vertex Radius " + tit, kTH2D, {axisPt, axisProdRadius}); - histos.add("particle/mothers/mothers/prodRadius3DVsPt", "Mothers mothers Prod. Vertex Radius XYZ " + tit, kTH2D, {axisPt, axisProdRadius}); + histos.add("particle/mothers/mothers/Number", "Number of Mothers mothers " + titlePdgCode, kTH1D, {{20, -0.5, 19.5}}); + histos.add("particle/mothers/mothers/PDGs", "Mothers mothers PDGs " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("particle/mothers/mothers/PDGsPrimaries", "Mothers mothers PDGs Primaries of " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("particle/mothers/mothers/PDGsSecondaries", "Mothers mothers PDGs Secondaries of " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("particle/mothers/mothers/prodVx", "Mothers mothers Prod. Vertex X " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/mothers/mothers/prodVy", "Mothers mothers Prod. Vertex Y " + titlePdgCode, kTH1D, {axes.axisProdXY}); + histos.add("particle/mothers/mothers/prodVz", "Mothers mothers Prod. Vertex Z " + titlePdgCode, kTH1D, {axes.axisProdZ}); + histos.add("particle/mothers/mothers/prodRadiusVsPt", "Mothers mothers Prod. Vertex Radius " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdRadius}); + histos.add("particle/mothers/mothers/prodRadius3DVsPt", "Mothers mothers Prod. Vertex Radius XYZ " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisProdRadius}); if (doprocessParticleOnly) { return; } histos.add("track/PDGs", "Track PDGs", kTH2D, {axisPDGs, axisCharge}); - histos.add("track/withoutParticle", "Tracks without particles", kTH1D, {axisPt}); + histos.add("track/withoutParticle", "Tracks without particles", kTH1D, {axes.axisMomentum}); histos.add("track/tofPDGs", "Track wTOF PDGs", kTH2D, {axisPDGs, axisCharge}); - histos.add("track/Pt", "Track Pt " + tit, kTH1D, {axisPt}); - histos.add("track/primariesPt", "Track Pt (primary) " + tit, kTH1D, {axisPt}); - histos.add("track/secondariesPt", "Track Pt (secondary) " + tit, kTH1D, {axisPt}); - histos.add("track/PtvsMCPt", "Track Pt vs MC Pt " + tit, kTH2D, {axisPt, axisPtMC}); - histos.add("track/Eta", "Track Eta " + tit, kTH1D, {axisEta}); - histos.add("track/primariesEta", "Track Eta (primary) " + tit, kTH1D, {axisEta}); - histos.add("track/secondariesEta", "Track Eta (secondary) " + tit, kTH1D, {axisEta}); - histos.add("track/primaries", "Source for primaries " + tit, kTH2D, {axisPDGs, axisCharge}); - histos.add("track/secondaries", "Source for secondaries " + tit, kTH2D, {axisPDGs, axisCharge}); + histos.add("track/Pt", "Track Pt " + titlePdgCode, kTH1D, {axes.axisMomentum}); + histos.add("track/primariesPt", "Track Pt (primary) " + titlePdgCode, kTH1D, {axes.axisMomentum}); + histos.add("track/secondariesPt", "Track Pt (secondary) " + titlePdgCode, kTH1D, {axes.axisMomentum}); + histos.add("track/PtvsMCPt", "Track Pt vs MC Pt " + titlePdgCode, kTH2D, {axes.axisMomentum, axes.axisMomentum}); + histos.add("track/Eta", "Track Eta " + titlePdgCode, kTH1D, {axes.axisEta}); + histos.add("track/primariesEta", "Track Eta (primary) " + titlePdgCode, kTH1D, {axes.axisEta}); + histos.add("track/secondariesEta", "Track Eta (secondary) " + titlePdgCode, kTH1D, {axes.axisEta}); + histos.add("track/primaries", "Source for primaries " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); + histos.add("track/secondaries", "Source for secondaries " + titlePdgCode, kTH2D, {axisPDGs, axisCharge}); } template @@ -230,236 +255,102 @@ struct Alice3SingleParticle { return Form("%i", particle.pdgCode()); } - void processStandard(const o2::aod::McCollisions& colls, - const soa::Join& tracks, - const aod::McParticles& mcParticles) + // Nch estimator: charged physical primaries in |eta| < settings.nchEtaMax + template + bool isNchParticle(ParticleType const& particle) { - for (const auto& col : colls) { - histos.fill(HIST("event/VtxX"), col.posX()); - histos.fill(HIST("event/VtxY"), col.posY()); - histos.fill(HIST("event/VtxZ"), col.posZ()); + if (!particle.isPhysicalPrimary() || std::abs(particle.eta()) >= settings.nchEtaMax) { + return false; } - std::vector ParticlesOfInterest; - for (const auto& mcParticle : mcParticles) { - const auto& pdgString = getPdgCodeString(mcParticle); - const auto& pdgCharge = getCharge(mcParticle); - histos.get(HIST("particle/PDGs"))->Fill(pdgString, pdgCharge, 1.f); - if (mcParticle.isPhysicalPrimary()) { - histos.get(HIST("particle/PDGsPrimaries"))->Fill(pdgString, pdgCharge, 1.f); - } else { - histos.get(HIST("particle/PDGsSecondaries"))->Fill(pdgString, pdgCharge, 1.f); - } - if (PDG != 0 && mcParticle.pdgCode() != PDG) { - continue; - } - if (mcParticle.y() < yMin || mcParticle.y() > yMax) { - continue; - } - if (mcParticle.eta() < etaMin || mcParticle.eta() > etaMax) { - continue; - } - if (mcParticle.pt() < ptMin || mcParticle.pt() > ptMax) { - continue; - } - if (mcParticle.vz() < prodMinZ || mcParticle.vz() > prodMaxZ) { - continue; - } - histos.fill(HIST("particle/producedProcess"), mcParticle.getProcess()); - histos.fill(HIST("particle/Pt"), mcParticle.pt()); - histos.fill(HIST("particle/P"), mcParticle.p()); - histos.fill(HIST("particle/Eta"), mcParticle.eta()); - if (mcParticle.isPhysicalPrimary()) { - histos.fill(HIST("particle/primariesPt"), mcParticle.pt()); - histos.fill(HIST("particle/primariesP"), mcParticle.p()); - histos.fill(HIST("particle/primariesEta"), mcParticle.eta()); - histos.fill(HIST("particle/primariesY"), mcParticle.y()); - } else { - histos.fill(HIST("particle/secondariesPt"), mcParticle.pt()); - histos.fill(HIST("particle/secondariesP"), mcParticle.p()); - histos.fill(HIST("particle/secondariesEta"), mcParticle.eta()); - histos.fill(HIST("particle/secondariesY"), mcParticle.y()); - } - histos.fill(HIST("particle/EvsPz"), mcParticle.e(), mcParticle.pz()); - histos.fill(HIST("particle/Y"), mcParticle.y()); - histos.fill(HIST("particle/YvzPz"), mcParticle.y(), mcParticle.pz()); - histos.fill(HIST("particle/EtavzPz"), mcParticle.eta(), mcParticle.pz()); - histos.fill(HIST("particle/PvzPz"), mcParticle.p(), mcParticle.pz()); - histos.fill(HIST("particle/PtvzPz"), mcParticle.pt(), mcParticle.pz()); - histos.fill(HIST("particle/Px"), mcParticle.px()); - histos.fill(HIST("particle/Py"), mcParticle.py()); - histos.fill(HIST("particle/Pz"), mcParticle.pz()); - histos.fill(HIST("particle/prodVx"), mcParticle.vx()); - histos.fill(HIST("particle/prodVy"), mcParticle.vy()); - histos.fill(HIST("particle/prodVz"), mcParticle.vz()); - if (mcParticle.has_daughters()) { - const auto& daughters = mcParticle.daughters_as(); - histos.fill(HIST("particle/daughters/Number"), daughters.size()); - bool filledDistance = false; - for (const auto& daughter : daughters) { - const auto& pdgStringDau = getPdgCodeString(daughter); - const auto& pdgChargeDau = getCharge(daughter); - - histos.get(HIST("particle/daughters/PDGs"))->Fill(pdgStringDau, pdgChargeDau, 1.f); - if (mcParticle.isPhysicalPrimary()) { - histos.get(HIST("particle/daughters/PDGsPrimaries"))->Fill(pdgStringDau, pdgChargeDau, 1.f); - } else { - histos.get(HIST("particle/daughters/PDGsSecondaries"))->Fill(pdgStringDau, pdgChargeDau, 1.f); - } - - histos.fill(HIST("particle/daughters/prodVx"), daughter.vx()); - histos.fill(HIST("particle/daughters/prodVy"), daughter.vy()); - histos.fill(HIST("particle/daughters/prodVz"), daughter.vz()); - if (filledDistance) { - continue; - } - filledDistance = true; - histos.fill(HIST("particle/daughters/prodRadiusVsPt"), mcParticle.pt(), std::sqrt(daughter.vx() * daughter.vx() + daughter.vy() * daughter.vy())); - histos.fill(HIST("particle/daughters/prodRadius3DVsPt"), mcParticle.pt(), std::sqrt(daughter.vx() * daughter.vx() + daughter.vy() * daughter.vy() + daughter.vz() * daughter.vz())); - - const float dx = daughter.vx() - mcParticle.vx(); - const float dy = daughter.vy() - mcParticle.vy(); - const float dz = daughter.vz() - mcParticle.vz(); - histos.fill(HIST("particle/decayDistanceX"), dx); - histos.fill(HIST("particle/decayDistanceY"), dy); - histos.fill(HIST("particle/decayDistanceZ"), dz); - histos.fill(HIST("particle/decayDistanceVsPt"), mcParticle.pt(), std::sqrt((dx) * (dx) + (dy) * (dy) + (dz) * (dz))); - if (pdg->GetParticle(mcParticle.pdgCode())) { - const float mass = pdg->GetParticle(mcParticle.pdgCode())->Mass(); - histos.fill(HIST("particle/decayCtau"), std::sqrt((dx) * (dx) + (dy) * (dy) + (dz) * (dz)) / mcParticle.pt() * mass * 10 * 1000); - } - } - } else { - histos.fill(HIST("particle/daughters/Number"), 0.f); - } - if (mcParticle.has_mothers()) { - const auto& mothers = mcParticle.mothers_as(); - for (const auto& mother : mothers) { - const auto& pdgStringMot = getPdgCodeString(mother); - const auto& pdgChargeMot = getCharge(mother); - - histos.get(HIST("particle/mothers/PDGs"))->Fill(pdgStringMot, pdgChargeMot, 1.f); - if (mcParticle.isPhysicalPrimary()) { - histos.get(HIST("particle/mothers/PDGsPrimaries"))->Fill(pdgStringMot, pdgChargeMot, 1.f); - } else { - histos.get(HIST("particle/mothers/PDGsSecondaries"))->Fill(pdgStringMot, pdgChargeMot, 1.f); - } - - histos.fill(HIST("particle/mothers/prodVx"), mother.vx()); - histos.fill(HIST("particle/mothers/prodVy"), mother.vy()); - histos.fill(HIST("particle/mothers/prodVz"), mother.vz()); - histos.fill(HIST("particle/mothers/prodRadiusVsPt"), mother.pt(), std::sqrt(mother.vx() * mother.vx() + mother.vy() * mother.vy())); - histos.fill(HIST("particle/mothers/prodRadius3DVsPt"), mother.pt(), std::sqrt(mother.vx() * mother.vx() + mother.vy() * mother.vy() + mother.vz() * mother.vz())); - if (mother.has_mothers()) { - const auto& mothers2 = mother.mothers_as(); - for (const auto& mother2 : mothers2) { - const auto& pdgStringMot2 = getPdgCodeString(mother2); - const auto& pdgChargeMot2 = getCharge(mother2); - - histos.get(HIST("particle/mothers/mothers/PDGs"))->Fill(pdgStringMot2, pdgChargeMot2, 1.f); - if (mcParticle.isPhysicalPrimary()) { - histos.get(HIST("particle/mothers/mothers/PDGsPrimaries"))->Fill(pdgStringMot2, pdgChargeMot2, 1.f); - } else { - histos.get(HIST("particle/mothers/mothers/PDGsSecondaries"))->Fill(pdgStringMot2, pdgChargeMot2, 1.f); - } - - histos.fill(HIST("particle/mothers/mothers/prodVx"), mother2.vx()); - histos.fill(HIST("particle/mothers/mothers/prodVy"), mother2.vy()); - histos.fill(HIST("particle/mothers/mothers/prodVz"), mother2.vz()); - histos.fill(HIST("particle/mothers/mothers/prodRadiusVsPt"), mother2.pt(), std::sqrt(mother2.vx() * mother2.vx() + mother2.vy() * mother2.vy())); - histos.fill(HIST("particle/mothers/mothers/prodRadius3DVsPt"), mother2.pt(), std::sqrt(mother2.vx() * mother2.vx() + mother2.vy() * mother2.vy() + mother2.vz() * mother2.vz())); - } - } - } - } + const auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); + return pdgParticle && pdgParticle->Charge() != 0; + } - histos.fill(HIST("particle/prodRadius"), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); - histos.fill(HIST("particle/prodVxVy"), mcParticle.vx(), mcParticle.vy()); - if (!mcParticle.isPhysicalPrimary()) { - if (mcParticle.getProcess() == 4) { - histos.fill(HIST("particle/prodVxVyStr"), mcParticle.vx(), mcParticle.vy()); - } else { - histos.fill(HIST("particle/prodVxVyMat"), mcParticle.vx(), mcParticle.vy()); - } - } else { - histos.fill(HIST("particle/prodVxVyPrm"), mcParticle.vx(), mcParticle.vy()); - } - histos.fill(HIST("particle/prodVxVsPt"), mcParticle.pt(), mcParticle.vx()); - histos.fill(HIST("particle/prodVyVsPt"), mcParticle.pt(), mcParticle.vy()); - histos.fill(HIST("particle/prodVzVsPt"), mcParticle.pt(), mcParticle.vz()); - histos.fill(HIST("particle/prodRadiusVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); - histos.fill(HIST("particle/prodRadiusVsEta"), mcParticle.eta(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); - histos.fill(HIST("particle/prodRadius3DVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy() + mcParticle.vz() * mcParticle.vz())); - ParticlesOfInterest.push_back(mcParticle.globalIndex()); + void fillSignalMultiplicity(const std::vector& nchPerCollision, const std::vector& nSignalPerCollision) + { + for (std::size_t iCollision = 0; iCollision < nchPerCollision.size(); ++iCollision) { + histos.fill(HIST("event/Nch"), nchPerCollision[iCollision]); + histos.fill(HIST("particle/NSignalVsNch"), nchPerCollision[iCollision], nSignalPerCollision[iCollision]); } + } + // Track histograms. Must be called after genericProcessParticles, which fills particlesOfInterest. + template + void genericProcessTracks(const TTracks& tracks) + { for (const auto& track : tracks) { if (!track.has_mcParticle()) { histos.fill(HIST("track/withoutParticle"), track.pt()); continue; } - const auto mcParticle = track.mcParticle_as(); + const auto mcParticle = track.template mcParticle_as(); histos.get(HIST("track/PDGs"))->Fill(getPdgCodeString(mcParticle), getCharge(mcParticle), 1.f); if (track.hasTOF()) { histos.get(HIST("track/tofPDGs"))->Fill(getPdgCodeString(mcParticle), getCharge(mcParticle), 1.f); } - if (IsStable.value == 0) { - LOG(info) << mcParticle.pdgCode() << " asked for " << PDG.value; + if (settings.isStable.value == 0) { + if (settings.doPrint) { + LOG(info) << mcParticle.pdgCode() << " asked for " << pdgCode.value; + } if (!mcParticle.has_mothers()) { continue; } // auto mothers = mcParticle.mothers(); - auto mothers = mcParticle.mothers_as(); - const auto ParticleIsInteresting = std::find(ParticlesOfInterest.begin(), ParticlesOfInterest.end(), mothers[0].globalIndex()) != ParticlesOfInterest.end(); - if (!ParticleIsInteresting) { + auto mothers = mcParticle.template mothers_as(); + const auto particleIsInteresting = std::find(particlesOfInterest.begin(), particlesOfInterest.end(), mothers[0].globalIndex()) != particlesOfInterest.end(); + if (!particleIsInteresting) { continue; } - if (doPrint) { + if (settings.doPrint) { LOG(info) << "Track " << track.globalIndex() << " comes from a " << mothers[0].pdgCode() << " and is a " << mcParticle.pdgCode(); } } else { - if (PDG != 0 && mcParticle.pdgCode() != PDG) { + if (pdgCode != 0 && mcParticle.pdgCode() != pdgCode) { continue; } - histos.fill(HIST("track/Pt"), track.pt() * charge); - histos.fill(HIST("track/PtvsMCPt"), track.pt() * charge, mcParticle.pt()); + histos.fill(HIST("track/Pt"), track.pt() * settings.charge); + histos.fill(HIST("track/PtvsMCPt"), track.pt() * settings.charge, mcParticle.pt()); histos.fill(HIST("track/Eta"), track.eta()); if (mcParticle.isPhysicalPrimary()) { - histos.fill(HIST("track/primariesPt"), track.pt() * charge); + histos.fill(HIST("track/primariesPt"), track.pt() * settings.charge); histos.fill(HIST("track/primariesEta"), track.eta()); } else { - histos.fill(HIST("track/secondariesPt"), track.pt() * charge); + histos.fill(HIST("track/secondariesPt"), track.pt() * settings.charge); histos.fill(HIST("track/secondariesEta"), track.eta()); } if (!mcParticle.has_mothers()) { - if (doPrint) { + if (settings.doPrint) { LOG(info) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode(); } continue; } // auto mothers = mcParticle.mothers(); - auto mothers = mcParticle.mothers_as(); + auto mothers = mcParticle.template mothers_as(); if (mcParticle.isPhysicalPrimary()) { histos.get(HIST("track/primaries"))->Fill(getPdgCodeString(mothers[0]), getCharge(mothers[0]), 1.f); } else { histos.get(HIST("track/secondaries"))->Fill(getPdgCodeString(mothers[0]), getCharge(mothers[0]), 1.f); } - if (doPrint) { + if (settings.doPrint) { LOG(info) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode() << " and comes from a " << mothers[0].pdgCode() << " and is " << (mcParticle.isPhysicalPrimary() ? "" : "not") << " a primary"; } } } } - PROCESS_SWITCH(Alice3SingleParticle, processStandard, "Process IU tracks", true); - void processParticleOnly(const o2::aod::McCollisions& colls, - const aod::McParticles& mcParticles) + // Event and MC particle histograms. Stores the selected MC particles in particlesOfInterest. + void genericProcessParticles(const o2::aod::McCollisions& colls, + const aod::McParticles& mcParticles) { for (const auto& col : colls) { histos.fill(HIST("event/VtxX"), col.posX()); histos.fill(HIST("event/VtxY"), col.posY()); histos.fill(HIST("event/VtxZ"), col.posZ()); } + std::vector nchPerCollision(colls.size(), 0); + std::vector nSignalPerCollision(colls.size(), 0); + particlesOfInterest.clear(); for (const auto& mcParticle : mcParticles) { const auto& pdgString = getPdgCodeString(mcParticle); const auto& pdgCharge = getCharge(mcParticle); @@ -469,18 +360,25 @@ struct Alice3SingleParticle { } else { histos.get(HIST("particle/PDGsSecondaries"))->Fill(pdgString, pdgCharge, 1.f); } - if (PDG != 0 && mcParticle.pdgCode() != PDG) { + if (isNchParticle(mcParticle)) { + ++nchPerCollision.at(mcParticle.mcCollisionId()); + } + if (pdgCode != 0 && mcParticle.pdgCode() != pdgCode) { continue; } - if (mcParticle.y() < yMin || mcParticle.y() > yMax) { + if (!yRange.contains(mcParticle.y())) { continue; } - if (mcParticle.eta() < etaMin || mcParticle.eta() > etaMax) { + if (!etaRange.contains(mcParticle.eta())) { continue; } - if (mcParticle.vz() < prodMinZ || mcParticle.vz() > prodMaxZ) { + if (!ptRange.contains(mcParticle.pt())) { continue; } + if (!prodZRange.contains(mcParticle.vz())) { + continue; + } + ++nSignalPerCollision.at(mcParticle.mcCollisionId()); histos.fill(HIST("particle/producedProcess"), mcParticle.getProcess()); histos.fill(HIST("particle/Pt"), mcParticle.pt()); histos.fill(HIST("particle/P"), mcParticle.p()); @@ -549,7 +447,7 @@ struct Alice3SingleParticle { histos.fill(HIST("particle/daughters/Number"), 0.f); } if (mcParticle.has_mothers()) { - auto mothers = mcParticle.mothers_as(); + const auto& mothers = mcParticle.mothers_as(); for (const auto& mother : mothers) { const auto& pdgStringMot = getPdgCodeString(mother); const auto& pdgChargeMot = getCharge(mother); @@ -566,13 +464,33 @@ struct Alice3SingleParticle { histos.fill(HIST("particle/mothers/prodVz"), mother.vz()); histos.fill(HIST("particle/mothers/prodRadiusVsPt"), mother.pt(), std::sqrt(mother.vx() * mother.vx() + mother.vy() * mother.vy())); histos.fill(HIST("particle/mothers/prodRadius3DVsPt"), mother.pt(), std::sqrt(mother.vx() * mother.vx() + mother.vy() * mother.vy() + mother.vz() * mother.vz())); + if (mother.has_mothers()) { + const auto& mothers2 = mother.mothers_as(); + for (const auto& mother2 : mothers2) { + const auto& pdgStringMot2 = getPdgCodeString(mother2); + const auto& pdgChargeMot2 = getCharge(mother2); + + histos.get(HIST("particle/mothers/mothers/PDGs"))->Fill(pdgStringMot2, pdgChargeMot2, 1.f); + if (mcParticle.isPhysicalPrimary()) { + histos.get(HIST("particle/mothers/mothers/PDGsPrimaries"))->Fill(pdgStringMot2, pdgChargeMot2, 1.f); + } else { + histos.get(HIST("particle/mothers/mothers/PDGsSecondaries"))->Fill(pdgStringMot2, pdgChargeMot2, 1.f); + } + + histos.fill(HIST("particle/mothers/mothers/prodVx"), mother2.vx()); + histos.fill(HIST("particle/mothers/mothers/prodVy"), mother2.vy()); + histos.fill(HIST("particle/mothers/mothers/prodVz"), mother2.vz()); + histos.fill(HIST("particle/mothers/mothers/prodRadiusVsPt"), mother2.pt(), std::sqrt(mother2.vx() * mother2.vx() + mother2.vy() * mother2.vy())); + histos.fill(HIST("particle/mothers/mothers/prodRadius3DVsPt"), mother2.pt(), std::sqrt(mother2.vx() * mother2.vx() + mother2.vy() * mother2.vy() + mother2.vz() * mother2.vz())); + } + } } } histos.fill(HIST("particle/prodRadius"), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); histos.fill(HIST("particle/prodVxVy"), mcParticle.vx(), mcParticle.vy()); if (!mcParticle.isPhysicalPrimary()) { - if (mcParticle.getProcess() == 4) { + if (mcParticle.getProcess() == TMCProcess::kPDecay) { histos.fill(HIST("particle/prodVxVyStr"), mcParticle.vx(), mcParticle.vy()); } else { histos.fill(HIST("particle/prodVxVyMat"), mcParticle.vx(), mcParticle.vy()); @@ -586,187 +504,35 @@ struct Alice3SingleParticle { histos.fill(HIST("particle/prodRadiusVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); histos.fill(HIST("particle/prodRadiusVsEta"), mcParticle.eta(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); histos.fill(HIST("particle/prodRadius3DVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy() + mcParticle.vz() * mcParticle.vz())); + particlesOfInterest.push_back(mcParticle.globalIndex()); } + fillSignalMultiplicity(nchPerCollision, nSignalPerCollision); + } + + void processStandard(const o2::aod::McCollisions& colls, + const soa::Join& tracks, + const aod::McParticles& mcParticles) + { + genericProcessParticles(colls, mcParticles); + genericProcessTracks(tracks); + } + + void processParticleOnly(const o2::aod::McCollisions& colls, + const aod::McParticles& mcParticles) + { + genericProcessParticles(colls, mcParticles); } - PROCESS_SWITCH(Alice3SingleParticle, processParticleOnly, "Process Particle only", false); void processNonIU(const o2::aod::McCollisions& colls, const soa::Join& tracks, const aod::McParticles& mcParticles) { - for (const auto& col : colls) { - histos.fill(HIST("event/VtxX"), col.posX()); - histos.fill(HIST("event/VtxY"), col.posY()); - histos.fill(HIST("event/VtxZ"), col.posZ()); - } - std::vector ParticlesOfInterest; - for (const auto& mcParticle : mcParticles) { - const auto& pdgString = getPdgCodeString(mcParticle); - const auto& pdgCharge = getCharge(mcParticle); - histos.get(HIST("particle/PDGs"))->Fill(pdgString, pdgCharge, 1.f); - if (mcParticle.isPhysicalPrimary()) { - histos.get(HIST("particle/PDGsPrimaries"))->Fill(pdgString, pdgCharge, 1.f); - } else { - histos.get(HIST("particle/PDGsSecondaries"))->Fill(pdgString, pdgCharge, 1.f); - } - if (PDG != 0 && mcParticle.pdgCode() != PDG) { - continue; - } - if (mcParticle.y() < yMin || mcParticle.y() > yMax) { - continue; - } - if (mcParticle.eta() < etaMin || mcParticle.eta() > etaMax) { - continue; - } - if (mcParticle.pt() < ptMin || mcParticle.pt() > ptMax) { - continue; - } - histos.fill(HIST("particle/Pt"), mcParticle.pt()); - histos.fill(HIST("particle/P"), mcParticle.p()); - histos.fill(HIST("particle/Eta"), mcParticle.eta()); - if (mcParticle.isPhysicalPrimary()) { - histos.fill(HIST("particle/primariesPt"), mcParticle.pt()); - histos.fill(HIST("particle/primariesP"), mcParticle.p()); - histos.fill(HIST("particle/primariesEta"), mcParticle.eta()); - histos.fill(HIST("particle/primariesY"), mcParticle.y()); - } else { - histos.fill(HIST("particle/secondariesPt"), mcParticle.pt()); - histos.fill(HIST("particle/secondariesP"), mcParticle.p()); - histos.fill(HIST("particle/secondariesEta"), mcParticle.eta()); - histos.fill(HIST("particle/secondariesY"), mcParticle.y()); - } - histos.fill(HIST("particle/EvsPz"), mcParticle.e(), mcParticle.pz()); - histos.fill(HIST("particle/Y"), mcParticle.y()); - histos.fill(HIST("particle/YvzPz"), mcParticle.y(), mcParticle.pz()); - histos.fill(HIST("particle/EtavzPz"), mcParticle.eta(), mcParticle.pz()); - histos.fill(HIST("particle/PvzPz"), mcParticle.p(), mcParticle.pz()); - histos.fill(HIST("particle/PtvzPz"), mcParticle.pt(), mcParticle.pz()); - histos.fill(HIST("particle/Px"), mcParticle.px()); - histos.fill(HIST("particle/Py"), mcParticle.py()); - histos.fill(HIST("particle/Pz"), mcParticle.pz()); - histos.fill(HIST("particle/prodVx"), mcParticle.vx()); - histos.fill(HIST("particle/prodVy"), mcParticle.vy()); - histos.fill(HIST("particle/prodVz"), mcParticle.vz()); - if (mcParticle.has_daughters()) { - auto daughters = mcParticle.daughters_as(); - histos.fill(HIST("particle/daughters/Number"), daughters.size()); - for (const auto& daughter : daughters) { - const auto& pdgStringDau = getPdgCodeString(daughter); - const auto& pdgChargeDau = getCharge(daughter); - - histos.get(HIST("particle/daughters/PDGs"))->Fill(pdgStringDau, pdgChargeDau, 1.f); - if (mcParticle.isPhysicalPrimary()) { - histos.get(HIST("particle/daughters/PDGsPrimaries"))->Fill(pdgStringDau, pdgChargeDau, 1.f); - } else { - histos.get(HIST("particle/daughters/PDGsSecondaries"))->Fill(pdgStringDau, pdgChargeDau, 1.f); - } - - histos.fill(HIST("particle/daughters/prodVx"), daughter.vx()); - histos.fill(HIST("particle/daughters/prodVy"), daughter.vy()); - histos.fill(HIST("particle/daughters/prodVz"), daughter.vz()); - histos.fill(HIST("particle/daughters/prodRadiusVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); - histos.fill(HIST("particle/daughters/prodRadius3DVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy() + mcParticle.vz() * mcParticle.vz())); - } - } else { - histos.fill(HIST("particle/daughters/Number"), 0.f); - } - if (mcParticle.has_mothers()) { - auto mothers = mcParticle.mothers_as(); - for (const auto& mother : mothers) { - const auto& pdgStringMot = getPdgCodeString(mother); - const auto& pdgChargeMot = getCharge(mother); - - histos.get(HIST("particle/mothers/PDGs"))->Fill(pdgStringMot, pdgChargeMot, 1.f); - if (mcParticle.isPhysicalPrimary()) { - histos.get(HIST("particle/mothers/PDGsPrimaries"))->Fill(pdgStringMot, pdgChargeMot, 1.f); - } else { - histos.get(HIST("particle/mothers/PDGsSecondaries"))->Fill(pdgStringMot, pdgChargeMot, 1.f); - } - - histos.fill(HIST("particle/mothers/prodVx"), mother.vx()); - histos.fill(HIST("particle/mothers/prodVy"), mother.vy()); - histos.fill(HIST("particle/mothers/prodVz"), mother.vz()); - histos.fill(HIST("particle/mothers/prodRadiusVsPt"), mother.pt(), std::sqrt(mother.vx() * mother.vx() + mother.vy() * mother.vy())); - histos.fill(HIST("particle/mothers/prodRadius3DVsPt"), mother.pt(), std::sqrt(mother.vx() * mother.vx() + mother.vy() * mother.vy() + mother.vz() * mother.vz())); - } - } - - histos.fill(HIST("particle/prodRadius"), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); - histos.fill(HIST("particle/prodVxVy"), mcParticle.vx(), mcParticle.vy()); - if (!mcParticle.isPhysicalPrimary()) { - if (mcParticle.getProcess() == 4) { - histos.fill(HIST("particle/prodVxVyStr"), mcParticle.vx(), mcParticle.vy()); - } else { - histos.fill(HIST("particle/prodVxVyMat"), mcParticle.vx(), mcParticle.vy()); - } - } else { - histos.fill(HIST("particle/prodVxVyPrm"), mcParticle.vx(), mcParticle.vy()); - } - histos.fill(HIST("particle/prodVxVsPt"), mcParticle.pt(), mcParticle.vx()); - histos.fill(HIST("particle/prodVyVsPt"), mcParticle.pt(), mcParticle.vy()); - histos.fill(HIST("particle/prodVzVsPt"), mcParticle.pt(), mcParticle.vz()); - histos.fill(HIST("particle/prodRadiusVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); - histos.fill(HIST("particle/prodRadiusVsEta"), mcParticle.eta(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy())); - histos.fill(HIST("particle/prodRadius3DVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy() + mcParticle.vz() * mcParticle.vz())); - ParticlesOfInterest.push_back(mcParticle.globalIndex()); - } - - for (const auto& track : tracks) { - if (!track.has_mcParticle()) { - histos.fill(HIST("track/withoutParticle"), track.pt()); - continue; - } - const auto mcParticle = track.mcParticle_as(); - histos.get(HIST("track/PDGs"))->Fill(getPdgCodeString(mcParticle), getCharge(mcParticle), 1.f); - if (track.hasTOF()) { - histos.get(HIST("track/tofPDGs"))->Fill(getPdgCodeString(mcParticle), getCharge(mcParticle), 1.f); - } - if (IsStable.value == 0) { - if (!mcParticle.has_mothers()) { - continue; - } - const auto& mothers = mcParticle.mothers_as(); - const auto& ParticleIsInteresting = std::find(ParticlesOfInterest.begin(), ParticlesOfInterest.end(), mothers[0].globalIndex()) != ParticlesOfInterest.end(); - if (!ParticleIsInteresting) { - continue; - } - if (doPrint) { - LOG(info) << "Track " << track.globalIndex() << " comes from a " << mothers[0].pdgCode() << " and is a " << mcParticle.pdgCode(); - } - } else { - if (PDG != 0 && mcParticle.pdgCode() != PDG) { - continue; - } - histos.fill(HIST("track/Pt"), track.pt() * charge); - histos.fill(HIST("track/PtvsMCPt"), track.pt() * charge, mcParticle.pt()); - histos.fill(HIST("track/Eta"), track.eta()); - if (mcParticle.isPhysicalPrimary()) { - histos.fill(HIST("track/primariesPt"), track.pt() * charge); - histos.fill(HIST("track/primariesEta"), track.eta()); - } else { - histos.fill(HIST("track/secondariesPt"), track.pt() * charge); - histos.fill(HIST("track/secondariesEta"), track.eta()); - } - - if (!mcParticle.has_mothers()) { - if (doPrint) { - LOG(info) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode(); - } - continue; - } - // auto mothers = mcParticle.mothers(); - auto mothers = mcParticle.mothers_as(); - if (mcParticle.isPhysicalPrimary()) { - histos.get(HIST("track/primaries"))->Fill(getPdgCodeString(mothers[0]), getCharge(mothers[0]), 1.f); - } else { - histos.get(HIST("track/secondaries"))->Fill(getPdgCodeString(mothers[0]), getCharge(mothers[0]), 1.f); - } - if (doPrint) { - LOG(info) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode() << " and comes from a " << mothers[0].pdgCode() << " and is " << (mcParticle.isPhysicalPrimary() ? "" : "not") << " a primary"; - } - } - } + genericProcessParticles(colls, mcParticles); + genericProcessTracks(tracks); } + + PROCESS_SWITCH(Alice3SingleParticle, processStandard, "Process IU tracks", false); + PROCESS_SWITCH(Alice3SingleParticle, processParticleOnly, "Process Particle only", true); PROCESS_SWITCH(Alice3SingleParticle, processNonIU, "Process non IU tracks", false); };