diff --git a/ALICE3/Core/TrackUtilities.h b/ALICE3/Core/TrackUtilities.h index 5ddce03fa15..edf92224edb 100644 --- a/ALICE3/Core/TrackUtilities.h +++ b/ALICE3/Core/TrackUtilities.h @@ -54,7 +54,7 @@ struct OTFParticle { // Getters int pdgCode() const { return mPdgCode; } - int isAlive() const { return mIsAlive; } + bool isAlive() const { return mIsAlive; } float vx() const { return mVx; } float vy() const { return mVy; } float vz() const { return mVz; } diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx index f8bcb033654..b365cd2f1fe 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -73,7 +73,7 @@ struct OnTheFlyDecayer { o2::upgrade::Decayer decayer; Service pdgDB; - std::map> mDecayDaughters; + std::map> mDecayDaughters; Configurable seed{"seed", 0, "Set seed for particle decayer"}; Configurable magneticField{"magneticField", 20., "Magnetic field (kG)"}; @@ -84,9 +84,62 @@ struct OnTheFlyDecayer { HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "pt axis for QA histograms"}; + struct McParticleAlice3 { + McParticleAlice3() = default; + ~McParticleAlice3() = default; + McParticleAlice3(const McParticleAlice3& src) = default; + McParticleAlice3(int collisionId, + int pdgCode, + int statusCode, + int flags, + int mother0, + int mother1, + int daughter0, + int daughter1, + float weight, + float px, float py, float pz, float e, + float vx, float vy, float vz, float vt, + float phi, float eta, float pt, float p, float y, + bool isAlive, bool isPrimary) : collisionId(collisionId), + pdgCode(pdgCode), + statusCode(statusCode), + flags(flags), + mothersIds{mother0, mother1}, + daughtersIdSlice{daughter0, daughter1}, + weight(weight), + px(px), + py(py), + pz(pz), + e(e), + vx(vx), + vy(vy), + vz(vz), + vt(vt), + phi(phi), + eta(eta), + pt(pt), + p(p), + y(y), + isAlive(isAlive), + isPrimary(isPrimary) {} + int collisionId; + int pdgCode; + int statusCode; + int flags; + int mothersIds[2]; + int daughtersIdSlice[2]; + float weight; + float px, py, pz, e; + float vx, vy, vz, vt; + float phi, eta, pt, p, y; + bool isAlive; + bool isPrimary; + }; + std::vector mEnabledDecays; void init(o2::framework::InitContext&) { + LOG(info) << "Initializing on-the-fly-decayer."; decayer.setSeed(seed); decayer.setBField(magneticField); for (int i = 0; i < NumDecays; ++i) { @@ -115,11 +168,13 @@ struct OnTheFlyDecayer { return std::find(mEnabledDecays.begin(), mEnabledDecays.end(), pdgCode) != mEnabledDecays.end(); } + std::vector mcParticlesAlice3; void process(aod::McCollision const&, aod::McParticles const& mcParticles) { mDecayDaughters.clear(); + mcParticlesAlice3.clear(); u_int64_t nStoredDaughters = 0; - for (int64_t index{0}; index < mcParticles.size(); ++index) { + for (int index{0}; index < static_cast(mcParticles.size()); ++index) { const auto& particle = mcParticles.iteratorAt(index); std::vector decayDaughters; static constexpr int MaxNestedDecays = 10; @@ -129,13 +184,14 @@ struct OnTheFlyDecayer { o2::upgrade::convertMCParticleToO2Track(particle, o2track, pdgDB); decayDaughters = decayer.decayParticle(pdgDB, o2track, particle.pdgCode()); for (size_t idau{0}; idau < decayDaughters.size(); ++idau) { - o2::upgrade::OTFParticle& dau = decayDaughters[idau]; + o2::upgrade::OTFParticle dau = decayDaughters[idau]; o2::track::TrackParCov dauTrack; o2::upgrade::convertOTFParticleToO2Track(dau, dauTrack, pdgDB); if (canDecay(dau.pdgCode())) { dau.setIsAlive(false); std::vector cascadingDaughers = decayer.decayParticle(pdgDB, dauTrack, dau.pdgCode()); - for (const auto& daudau : cascadingDaughers) { + for (size_t idaudau{0}; idaudau < cascadingDaughers.size(); ++idaudau) { + o2::upgrade::OTFParticle daudau = cascadingDaughers[idaudau]; decayDaughters.push_back(daudau); if (MaxNestedDecays < ++nDecays) { LOG(error) << "Seemingly stuck trying to perpetually decay products from pdg: " << particle.pdgCode(); @@ -194,14 +250,13 @@ struct OnTheFlyDecayer { daughtersIdSlice[1] = static_cast(particle.daughtersIds()[1]); } - std::span motherSpan(particle.mothersIds().data(), particle.mothersIds().size()); mDecayDaughters.emplace(index, decayDaughters); nStoredDaughters += decayDaughters.size(); - float phi = o2::constants::math::PI + std::atan2(-1.0f * particle.py(), -1.0f * particle.px()); + const float phi = o2::constants::math::PI + std::atan2(-1.0f * particle.py(), -1.0f * particle.px()); float eta; // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1922 - float pt = std::sqrt(particle.px() * particle.px() + particle.py() * particle.py()); - float p = std::sqrt(particle.px() * particle.px() + particle.py() * particle.py() + particle.pz() * particle.pz()); + const float pt = std::sqrt(particle.px() * particle.px() + particle.py() * particle.py()); + const float p = std::sqrt(particle.px() * particle.px() + particle.py() * particle.py() + particle.pz() * particle.pz()); float y; // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943 if ((p - particle.pz()) < Tolerance) { @@ -218,29 +273,29 @@ struct OnTheFlyDecayer { // TODO: Particle status code // TODO: Expression columns - tableMcParticlesWithDau(particle.mcCollisionId(), particle.pdgCode(), particle.statusCode(), - particle.flags(), motherSpan, daughtersIdSlice, particle.weight(), - particle.px(), particle.py(), particle.pz(), particle.e(), - particle.vx(), particle.vy(), particle.vz(), particle.vt(), - phi, eta, pt, p, y, !canDecay(particle.pdgCode()), true); + auto mothers = particle.mothersIds(); + int mother0 = mothers.size() > 0 ? mothers[0] : -1; + int mother1 = mothers.size() > 1 ? mothers[1] : mother0; + mcParticlesAlice3.push_back(McParticleAlice3{particle.mcCollisionId(), particle.pdgCode(), particle.statusCode(), + particle.flags(), mother0, mother1, + daughtersIdSlice[0], daughtersIdSlice[1], particle.weight(), + particle.px(), particle.py(), particle.pz(), particle.e(), + particle.vx(), particle.vy(), particle.vz(), particle.vt(), + phi, eta, pt, p, y, !canDecay(particle.pdgCode()), true}); } int daughtersIdSlice[2] = {-1, -1}; for (const auto& [index, decayDaughters] : mDecayDaughters) { for (const auto& dau : decayDaughters) { if (index >= mcParticles.size()) { - LOG(warn) << "--- Index " << index << " out of bounds for mcParticles table of size " << mcParticles.size(); + LOG(error) << "--- Index " << index << " out of bounds for mcParticles table of size " << mcParticles.size() << std::endl; continue; } - auto mother = mcParticles.iteratorAt(index); - std::vector motherIds = {static_cast(index)}; - std::span motherSpan(motherIds.data(), motherIds.size()); - - float phi = o2::constants::math::PI + std::atan2(-1.0f * dau.py(), -1.0f * dau.px()); + const float phi = o2::constants::math::PI + std::atan2(-1.0f * dau.py(), -1.0f * dau.px()); float eta; // Conditional as https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1922 - float pt = std::sqrt(dau.px() * dau.px() + dau.py() * dau.py()); - float p = std::sqrt(dau.px() * dau.px() + dau.py() * dau.py() + dau.pz() * dau.pz()); + const float pt = std::sqrt(dau.px() * dau.px() + dau.py() * dau.py()); + const float p = std::sqrt(dau.px() * dau.px() + dau.py() * dau.py() + dau.pz() * dau.pz()); float y; // Conditional as https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943 if ((p - dau.pz()) < Tolerance) { @@ -283,13 +338,25 @@ struct OnTheFlyDecayer { // TODO: Particle status code // TODO: Expression columns // TODO: vt - tableMcParticlesWithDau(mother.mcCollisionId(), dau.pdgCode(), 1, - mother.flags(), motherSpan, daughtersIdSlice, mother.weight(), - dau.px(), dau.py(), dau.pz(), dau.e(), - dau.vx(), dau.vy(), dau.vz(), mother.vt(), - phi, eta, pt, p, y, dau.isAlive(), false); + auto mother = mcParticles.iteratorAt(index); + mcParticlesAlice3.push_back(McParticleAlice3{mother.mcCollisionId(), dau.pdgCode(), 1, + -1, index, index, daughtersIdSlice[0], daughtersIdSlice[1], mother.weight(), + dau.px(), dau.py(), dau.pz(), dau.e(), + dau.vx(), dau.vy(), dau.vz(), mother.vt(), + phi, eta, pt, p, y, dau.isAlive(), false}); } } + + for (const auto& particle : mcParticlesAlice3) { + std::span motherSpan(particle.mothersIds, 2); + + tableMcParticlesWithDau(particle.collisionId, particle.pdgCode, particle.statusCode, + particle.flags, motherSpan, particle.daughtersIdSlice, particle.weight, + particle.px, particle.py, particle.pz, particle.e, + particle.vx, particle.vy, particle.vz, particle.vt, + particle.phi, particle.eta, particle.pt, particle.p, particle.y, + particle.isAlive, particle.isPrimary); + } } }; diff --git a/ALICE3/Tasks/alice3DecayerQA.cxx b/ALICE3/Tasks/alice3DecayerQA.cxx index 7eff59fe2d0..5009df9f437 100644 --- a/ALICE3/Tasks/alice3DecayerQA.cxx +++ b/ALICE3/Tasks/alice3DecayerQA.cxx @@ -36,7 +36,23 @@ using namespace o2::framework; struct Alice3DecayerQA { HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "pt axis for QA histograms"}; + + struct : ConfigurableGroup { + ConfigurableAxis axisCollisionId{"axisCollisionId", {1000, 0, 999}, "CollisionId axis for QA histograms"}; + ConfigurableAxis axisPdgCode{"axisPdgCode", {1000, 0, 999}, "PdgCode axis for QA histograms"}; + ConfigurableAxis axisStatusCode{"axisStatusCode", {1000, 0, 999}, "StatusCode axis for QA histograms"}; + ConfigurableAxis axisFlags{"axisFlags", {10, 0, 9}, "Flags axis for QA histograms"}; + ConfigurableAxis axisMothersIds{"axisMothersIds", {1000, 0, 999}, "MothersIds axis for QA histograms"}; + ConfigurableAxis axisDaughtersIds{"axisDaughtersIds", {1000, 0, 999}, "DaughtersIds axis for QA histograms"}; + ConfigurableAxis axisWeight{"axisWeight", {2, 0, 1}, "Weight axis for QA histograms"}; + ConfigurableAxis axisPos{"axisPos", {1000, 0, 999}, "Position axis for QA histograms"}; + ConfigurableAxis axisPhi{"axisPhi", {720, -360, 360}, "Phi axis for QA histograms"}; + ConfigurableAxis axisEta{"axisEta", {80, -4, 4}, "Eta axis for QA histograms"}; + ConfigurableAxis axisRapidity{"axisRapidity", {80, -4, 4}, "Rapidity axis for QA histograms"}; + ConfigurableAxis axisIsAlive{"axisIsAlive", {2, 0, 1}, "IsAlive axis for QA histograms"}; + ConfigurableAxis axisIsPrimary{"axisIsPrimary", {2, 0, 1}, "IsPrimary axis for QA histograms"}; + ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "pt axis for QA histograms"}; + } axes; Partition trueEl = aod::mcparticle::pdgCode == static_cast(kElectron); Partition trueMu = aod::mcparticle::pdgCode == static_cast(kMuonMinus); @@ -52,58 +68,110 @@ struct Alice3DecayerQA { void init(o2::framework::InitContext&) { - histos.add("DefaultMC/hElPt", "hElPt", kTH1D, {axisPt}); - histos.add("DefaultMC/hMuPt", "hMuPt", kTH1D, {axisPt}); - histos.add("DefaultMC/hPiPt", "hPiPt", kTH1D, {axisPt}); - histos.add("DefaultMC/hKaPt", "hKaPt", kTH1D, {axisPt}); - histos.add("DefaultMC/hPrPt", "hPrPt", kTH1D, {axisPt}); - - histos.add("MCWithDau/hElPt", "hElPt", kTH1D, {axisPt}); - histos.add("MCWithDau/hMuPt", "hMuPt", kTH1D, {axisPt}); - histos.add("MCWithDau/hPiPt", "hPiPt", kTH1D, {axisPt}); - histos.add("MCWithDau/hKaPt", "hKaPt", kTH1D, {axisPt}); - histos.add("MCWithDau/hPrPt", "hPrPt", kTH1D, {axisPt}); + histos.add("DefaultMC/hElPt", "hElPt", kTH1D, {axes.axisPt}); + histos.add("DefaultMC/hMuPt", "hMuPt", kTH1D, {axes.axisPt}); + histos.add("DefaultMC/hPiPt", "hPiPt", kTH1D, {axes.axisPt}); + histos.add("DefaultMC/hKaPt", "hKaPt", kTH1D, {axes.axisPt}); + histos.add("DefaultMC/hPrPt", "hPrPt", kTH1D, {axes.axisPt}); + + histos.add("MCWithDau/hElPt", "hElPt", kTH1D, {axes.axisPt}); + histos.add("MCWithDau/hMuPt", "hMuPt", kTH1D, {axes.axisPt}); + histos.add("MCWithDau/hPiPt", "hPiPt", kTH1D, {axes.axisPt}); + histos.add("MCWithDau/hKaPt", "hKaPt", kTH1D, {axes.axisPt}); + histos.add("MCWithDau/hPrPt", "hPrPt", kTH1D, {axes.axisPt}); + + histos.add("MCWithDau/hCollisionId", "hCollisionId", kTH1D, {axes.axisCollisionId}); + histos.add("MCWithDau/hPdgCode", "hPdgCode", kTH1D, {axes.axisPdgCode}); + histos.add("MCWithDau/hStatusCode", "hStatusCode", kTH1D, {axes.axisStatusCode}); + histos.add("MCWithDau/hFlags", "hFlags", kTH1D, {axes.axisFlags}); + histos.add("MCWithDau/hMothersIds", "hMothersIds", kTH1D, {axes.axisMothersIds}); + histos.add("MCWithDau/hDaughtersIds", "hDaughtersIds", kTH1D, {axes.axisDaughtersIds}); + histos.add("MCWithDau/hWeight", "hWeight", kTH1D, {axes.axisWeight}); + histos.add("MCWithDau/hVx", "hVx", kTH1D, {axes.axisPos}); + histos.add("MCWithDau/hVy", "hVy", kTH1D, {axes.axisPos}); + histos.add("MCWithDau/hVz", "hVz", kTH1D, {axes.axisPos}); + histos.add("MCWithDau/hVt", "hVt", kTH1D, {axes.axisPos}); + histos.add("MCWithDau/hPhi", "hPhi", kTH1D, {axes.axisPhi}); + histos.add("MCWithDau/hEta", "hEta", kTH1D, {axes.axisEta}); + histos.add("MCWithDau/hRapidity", "hRapidity", kTH1D, {axes.axisRapidity}); + histos.add("MCWithDau/hIsAlive", "hIsAlive", kTH1D, {axes.axisIsAlive}); + histos.add("MCWithDau/hIsPrimary", "hIsPrimary", kTH1D, {axes.axisIsPrimary}); + histos.add("MCWithDau/hPx", "hPx", kTH1D, {axes.axisPt}); + histos.add("MCWithDau/hPy", "hPy", kTH1D, {axes.axisPt}); + histos.add("MCWithDau/hPz", "hPz", kTH1D, {axes.axisPt}); + histos.add("MCWithDau/hPt", "hPt", kTH1D, {axes.axisPt}); + histos.add("MCWithDau/hP", "hP", kTH1D, {axes.axisPt}); + histos.add("MCWithDau/hE", "hE", kTH1D, {axes.axisPt}); } void processMC(const aod::McParticles&) { - for (auto const& particle : trueEl) { + for (const auto& particle : trueEl) { histos.fill(HIST("DefaultMC/hElPt"), particle.pt()); } - for (auto const& particle : trueMu) { + for (const auto& particle : trueMu) { histos.fill(HIST("DefaultMC/hMuPt"), particle.pt()); } - for (auto const& particle : truePi) { + for (const auto& particle : truePi) { histos.fill(HIST("DefaultMC/hPiPt"), particle.pt()); } - for (auto const& particle : trueKa) { + for (const auto& particle : trueKa) { histos.fill(HIST("DefaultMC/hKaPt"), particle.pt()); } - for (auto const& particle : truePr) { + for (const auto& particle : truePr) { histos.fill(HIST("DefaultMC/hPrPt"), particle.pt()); } } - void processMCWithDau(const aod::McParticlesWithDau&) + void processMCWithDau(const aod::McCollision&, const aod::McParticlesWithDau& particles) { - for (auto const& particle : trueElWithDau) { + for (const auto& particle : trueElWithDau) { histos.fill(HIST("MCWithDau/hElPt"), particle.pt()); } - for (auto const& particle : trueMuWithDau) { + for (const auto& particle : trueMuWithDau) { histos.fill(HIST("MCWithDau/hMuPt"), particle.pt()); } - for (auto const& particle : truePiWithDau) { + for (const auto& particle : truePiWithDau) { histos.fill(HIST("MCWithDau/hPiPt"), particle.pt()); } - for (auto const& particle : trueKaWithDau) { + for (const auto& particle : trueKaWithDau) { histos.fill(HIST("MCWithDau/hKaPt"), particle.pt()); } - for (auto const& particle : truePrWithDau) { + for (const auto& particle : truePrWithDau) { histos.fill(HIST("MCWithDau/hPrPt"), particle.pt()); } + + for (const auto& particle : particles) { + histos.fill(HIST("MCWithDau/hCollisionId"), particle.mcCollisionId()); + histos.fill(HIST("MCWithDau/hPdgCode"), particle.pdgCode()); + histos.fill(HIST("MCWithDau/hStatusCode"), particle.statusCode()); + histos.fill(HIST("MCWithDau/hFlags"), particle.flags()); + histos.fill(HIST("MCWithDau/hWeight"), particle.weight()); + histos.fill(HIST("MCWithDau/hVx"), particle.vx()); + histos.fill(HIST("MCWithDau/hVy"), particle.vy()); + histos.fill(HIST("MCWithDau/hVz"), particle.vz()); + histos.fill(HIST("MCWithDau/hVt"), particle.vt()); + histos.fill(HIST("MCWithDau/hPhi"), particle.phi()); + histos.fill(HIST("MCWithDau/hEta"), particle.eta()); + histos.fill(HIST("MCWithDau/hRapidity"), particle.y()); + histos.fill(HIST("MCWithDau/hIsAlive"), particle.isAlive()); + histos.fill(HIST("MCWithDau/hIsPrimary"), particle.isPrimary()); + histos.fill(HIST("MCWithDau/hPx"), particle.px()); + histos.fill(HIST("MCWithDau/hPy"), particle.py()); + histos.fill(HIST("MCWithDau/hPz"), particle.pz()); + histos.fill(HIST("MCWithDau/hPt"), particle.pt()); + histos.fill(HIST("MCWithDau/hP"), particle.p()); + histos.fill(HIST("MCWithDau/hE"), particle.e()); + for (const auto& motherParticleId : particle.mothersIds()) { + histos.fill(HIST("MCWithDau/hMothersIds"), motherParticleId); + } + for (const auto& dauParticleId : particle.mothersIds()) { + histos.fill(HIST("MCWithDau/hDaughtersIds"), dauParticleId); + } + } } - PROCESS_SWITCH(Alice3DecayerQA, processMC, "fill MC-only histograms", true); + PROCESS_SWITCH(Alice3DecayerQA, processMC, "fill MC-only histograms", false); PROCESS_SWITCH(Alice3DecayerQA, processMCWithDau, "fill MC-with-dau histograms", true); };