From 8973025727e3d465f9827c9c1cad40642255fd75 Mon Sep 17 00:00:00 2001 From: jesgum Date: Thu, 18 Dec 2025 14:04:16 +0100 Subject: [PATCH 1/7] First step towards general handling of otf decays --- ALICE3/Core/Decayer.h | 140 +++++++++++++++++++ ALICE3/Core/TrackUtilities.h | 32 ++++- ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx | 129 +++++++++++++++++ 3 files changed, 300 insertions(+), 1 deletion(-) create mode 100644 ALICE3/Core/Decayer.h create mode 100644 ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h new file mode 100644 index 00000000000..9323a2c0e49 --- /dev/null +++ b/ALICE3/Core/Decayer.h @@ -0,0 +1,140 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file Decayer.h +/// \author Jesper Karlsson Gumprecht +/// \since 15/12/2025 +/// + +#ifndef ALICE3_CORE_DECAYER_H_ +#define ALICE3_CORE_DECAYER_H_ + +#include "ALICE3/Core/TrackUtilities.h" +#include "ReconstructionDataFormats/Track.h" + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +namespace o2 +{ +namespace upgrade +{ + +class Decayer { + public: + // Default constructor + Decayer() = default; + + template + std::vector decayParticle(const TDatabase& pdgDB, const o2::track::TrackParCov& track, const int pdgCode) + { + const auto& particleInfo = pdgDB->GetParticle(pdgCode); + const int charge = particleInfo->Charge() / 3; + const double mass = particleInfo->Mass(); + std::array mom; + std::array pos; + track.getPxPyPzGlo(mom); + track.getXYZGlo(pos); + + const double u = mRand3.Uniform(0, 1); + const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm + const double betaGamma = track.getP() / mass; + const double rxyz = -betaGamma * ctau * std::log(1 - u); + + float sna, csa; + o2::math_utils::CircleXYf_t circle; + track.getCircleParams(mBz, circle, sna, csa); + const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl()); + const double theta = rxy / circle.rC; + + const double vx = ((pos[0] - circle.xC) * std::cos(theta) - (pos[1] - circle.yC) * std::sin(theta)) + circle.xC; + const double vy = ((pos[1] - circle.yC) * std::cos(theta) + (pos[0] - circle.xC) * std::sin(theta)) + circle.yC; + const double vz = mom[2] + rxyz * (mom[2] / track.getP()); + + const double px = mom[0] * std::cos(theta) - mom[1] * std::sin(theta); + const double py = mom[1] * std::cos(theta) + mom[0] * std::sin(theta); + const double e = std::sqrt(mass * mass + px * px + py * py + mom[2] * mom[2]); + + double brTotal = 0.; + for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) { + brTotal += particleInfo->DecayChannel(ch)->BranchingRatio(); + } + + double brSum = 0.; + std::vector dauMasses; + std::vector pdgCodesDaughters; + const double randomChannel = mRand3.Uniform(0., brTotal); + for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) { + brSum += particleInfo->DecayChannel(ch)->BranchingRatio(); + if (randomChannel < brSum) { + for (int dau = 0; dau < particleInfo->DecayChannel(ch)->NDaughters(); ++dau) { + const int pdgDau = particleInfo->DecayChannel(ch)->DaughterPdgCode(dau); + pdgCodesDaughters.push_back(pdgDau); + const auto& dauInfo = pdgDB->GetParticle(pdgDau); + dauMasses.push_back(dauInfo->Mass()); + } + break; + } + } + + TLorentzVector tlv(px, py, mom[2], e); + TGenPhaseSpace decay; + decay.SetDecay(tlv, dauMasses.size(), dauMasses.data()); + decay.Generate(); + + std::vector decayProducts; + for (size_t i = 0; i < dauMasses.size(); ++i) { + o2::upgrade::OTFParticle particle; + TLorentzVector dau = *decay.GetDecay(i); + particle.setPDG(pdgCodesDaughters[i]); + particle.setVxVyVz(vx, vy, vz); + particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E()); + decayProducts.push_back(particle); + } + + return decayProducts; + } + + + // Setters + void setSeed(const int seed) + { + mRand3.SetSeed(seed); // For decay length sampling + gRandom->SetSeed(seed); // For TGenPhaseSpace + } + + void setBField(const double b) { mBz = b; } + + // Getters + + + + private: + TRandom3 mRand3; + double mBz; +}; + +} // namespace fastsim +} // namespace o2 + +#endif // ALICE3_CORE_DECAYER_H_ diff --git a/ALICE3/Core/TrackUtilities.h b/ALICE3/Core/TrackUtilities.h index 9bbc00c8f09..9543a71bc9a 100644 --- a/ALICE3/Core/TrackUtilities.h +++ b/ALICE3/Core/TrackUtilities.h @@ -21,9 +21,10 @@ #include "ReconstructionDataFormats/Track.h" #include "TLorentzVector.h" - #include + + namespace o2::upgrade { @@ -37,6 +38,17 @@ void convertTLorentzVectorToO2Track(const int charge, const std::vector productionVertex, o2::track::TrackParCov& o2track); +struct OTFParticle { + int pdgCode; + float e; + float vx, vy, vz; + float px, py, pz; + + void setPDG(int pdg) { pdgCode = pdg; } + void setVxVyVz(float _vx, float _vy, float _vz) { vx = _vx; vy = _vy; vz = _vz; } + void setPxPyPzE(float _px, float _py, float _pz, float _e) { px = _px; py = _py; pz = _pz; e = _e; } +}; + /// Function to convert a TLorentzVector into a perfect Track /// \param pdgCode particle pdg /// \param particle the particle to convert (TLorentzVector) @@ -58,6 +70,24 @@ void convertTLorentzVectorToO2Track(int pdgCode, convertTLorentzVectorToO2Track(charge, particle, productionVertex, o2track); } +/// Function to convert a OTFParticle into a perfect Track +/// \param particle the particle to convert (OTFParticle) +/// \param o2track the address of the resulting TrackParCov +/// \param pdg the pdg service +template +void convertOTFParticleToO2Track(const OTFParticle& particle, o2::track::TrackParCov& o2track, const PdgService& pdg) +{ + const auto pdgInfo = pdg->GetParticle(particle.pdgCode); + int charge = 0; + if (pdgInfo != nullptr) { + charge = pdgInfo->Charge() / 3; + } + + static TLorentzVector tlv; + tlv.SetPxPyPzE(particle.px, particle.py, particle.pz, particle.e); + convertTLorentzVectorToO2Track(particle.pdgCode, tlv, {particle.vx, particle.vy, particle.vz}, o2track, pdg); +} + /// Function to convert a McParticle into a perfect Track /// \param particle the particle to convert (mcParticle) /// \param o2track the address of the resulting TrackParCov diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx new file mode 100644 index 00000000000..b67bc2e7d92 --- /dev/null +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -0,0 +1,129 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file onTheFlyDecayer.cxx +/// \brief pre-processing for on-the-fly analysis +/// \author Jesper Karlsson Gumprecht +/// + +#include "ALICE3/Core/Decayer.h" +#include "ALICE3/Core/TrackUtilities.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +using namespace o2; +using namespace o2::framework; + +static constexpr int kNumDecays = 7; +static constexpr int kNumParameters = 1; +static constexpr int defaultParameters[kNumDecays][kNumParameters]{{1}, {1}, {1}, {1}, {1}, {1}, {1}}; +static const std::vector parameterNames{"enable"}; +static const std::vector particleNames{"K0s", + "Lambda", + "Anti-Lambda", + "Xi", + "Anti-Xi", + "Omega", + "Anti-Omega"}; +static const std::vector pdgCodes{kK0Short, + kLambda0, + kLambda0Bar, + kXiMinus, + kXiPlusBar, + kOmegaMinus, + kOmegaPlusBar}; + +struct OnTheFlyDecayer { + o2::upgrade::Decayer decayer; + Service ccdb; + Service pdgDB; + + Configurable seed{"seed", 0, "Set seed for particle decayer"}; + Configurable magneticField{"magneticField", 20., "Magnetic field (kG)"}; + Configurable> enabledDecays{"enabledDecays", + {defaultParameters[0], kNumDecays, kNumParameters, particleNames, parameterNames}, + "Enable option for particle to be decayed: 0 - no, 1 - yes"}; + + std::vector mEnabledDecays; + void init(o2::framework::InitContext&) + { + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setTimestamp(-1); + decayer.setSeed(seed); + decayer.setBField(magneticField); + for (int i = 0; i < kNumDecays; ++i) { + if (enabledDecays->get(particleNames[i].c_str(), "enable")) { + mEnabledDecays.push_back(pdgCodes[i]); + } + } + } + + bool canDecay(const int pdgCode) + { + return std::find(mEnabledDecays.begin(), mEnabledDecays.end(), pdgCode) != mEnabledDecays.end(); + } + + void process(aod::McCollision const&, aod::McParticles const& mcParticles) + { + for (const auto& particle : mcParticles) { + if (!canDecay(particle.pdgCode())) { + continue; + } + + o2::track::TrackParCov o2track; + o2::upgrade::convertMCParticleToO2Track(particle, o2track, pdgDB); + std::vector decayDaughters = decayer.decayParticle(pdgDB, o2track, particle.pdgCode()); + std::cout << particle.pdgCode() << " decayed into: "; + for (const auto& dau : decayDaughters) { + o2::track::TrackParCov dauTrack; + o2::upgrade::convertOTFParticleToO2Track(dau, dauTrack, pdgDB); + if (canDecay(dau.pdgCode)) { + std::vector cascadingDaughers = decayer.decayParticle(pdgDB, dauTrack, dau.pdgCode); + for (const auto& daudau : cascadingDaughers) { + decayDaughters.push_back(daudau); + } + } + } + std::vector isFinalState; + for (size_t i = 0; i < decayDaughters.size(); ++i) { + isFinalState.push_back(!canDecay(decayDaughters[i].pdgCode)); + } + } + } +}; + + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} From 42e84040d12ac5e9d5d9c948a60bb645ced946fb Mon Sep 17 00:00:00 2001 From: jesgum Date: Thu, 18 Dec 2025 21:33:51 +0100 Subject: [PATCH 2/7] Do v0 correctly --- ALICE3/Core/Decayer.h | 38 ++++++++------ ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx | 52 ++++++++++++++++++-- 2 files changed, 72 insertions(+), 18 deletions(-) diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h index 9323a2c0e49..efe6ecc768f 100644 --- a/ALICE3/Core/Decayer.h +++ b/ALICE3/Core/Decayer.h @@ -60,22 +60,32 @@ class Decayer { const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm const double betaGamma = track.getP() / mass; const double rxyz = -betaGamma * ctau * std::log(1 - u); - - float sna, csa; - o2::math_utils::CircleXYf_t circle; - track.getCircleParams(mBz, circle, sna, csa); - const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl()); - const double theta = rxy / circle.rC; - - const double vx = ((pos[0] - circle.xC) * std::cos(theta) - (pos[1] - circle.yC) * std::sin(theta)) + circle.xC; - const double vy = ((pos[1] - circle.yC) * std::cos(theta) + (pos[0] - circle.xC) * std::sin(theta)) + circle.yC; - const double vz = mom[2] + rxyz * (mom[2] / track.getP()); - - const double px = mom[0] * std::cos(theta) - mom[1] * std::sin(theta); - const double py = mom[1] * std::cos(theta) + mom[0] * std::sin(theta); - const double e = std::sqrt(mass * mass + px * px + py * py + mom[2] * mom[2]); + double vx, vy, vz; + double px, py, e; + + if (!charge) { + vx = pos[0] + rxyz * (mom[0] / track.getP()); + vy = pos[1] + rxyz * (mom[1] / track.getP()); + vy = pos[2] + rxyz * (mom[2] / track.getP()); + px = mom[0]; + py = mom[2]; + } else { + float sna, csa; + o2::math_utils::CircleXYf_t circle; + track.getCircleParams(mBz, circle, sna, csa); + const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl()); + const double theta = rxy / circle.rC; + + vx = ((pos[0] - circle.xC) * std::cos(theta) - (pos[1] - circle.yC) * std::sin(theta)) + circle.xC; + vy = ((pos[1] - circle.yC) * std::cos(theta) + (pos[0] - circle.xC) * std::sin(theta)) + circle.yC; + vz = mom[2] + rxyz * (mom[2] / track.getP()); + + px = mom[0] * std::cos(theta) - mom[1] * std::sin(theta); + py = mom[1] * std::cos(theta) + mom[0] * std::sin(theta); + } double brTotal = 0.; + e = std::sqrt(mass * mass + px * px + py * py + mom[2] * mom[2]); for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) { brTotal += particleInfo->DecayChannel(ch)->BranchingRatio(); } diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx index b67bc2e7d92..391e2096a82 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -17,6 +17,7 @@ #include "ALICE3/Core/Decayer.h" #include "ALICE3/Core/TrackUtilities.h" +#include "ALICE3/DataModel/OTFMCParticle.h" #include #include @@ -64,16 +65,23 @@ static const std::vector pdgCodes{kK0Short, kOmegaPlusBar}; struct OnTheFlyDecayer { + Produces tableMcParticles; o2::upgrade::Decayer decayer; Service ccdb; Service pdgDB; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; Configurable seed{"seed", 0, "Set seed for particle decayer"}; Configurable magneticField{"magneticField", 20., "Magnetic field (kG)"}; Configurable> enabledDecays{"enabledDecays", {defaultParameters[0], kNumDecays, kNumParameters, particleNames, parameterNames}, "Enable option for particle to be decayed: 0 - no, 1 - yes"}; - + + + ConfigurableAxis axisRadius{"axisRadius", {1000, 0, 100}, "Radius"}; + 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"}; + + std::vector mEnabledDecays; void init(o2::framework::InitContext&) { @@ -86,6 +94,14 @@ struct OnTheFlyDecayer { mEnabledDecays.push_back(pdgCodes[i]); } } + + histos.add("K0S/hGenK0S", "hGenK0S", kTH2D, {axisRadius, axisPt}); + histos.add("Lambda/hGenLambda", "hGenLambda", kTH2D, {axisRadius, axisPt}); + histos.add("AntiLambda/hGenAntiLambda", "hGenAntiLambda", kTH2D, {axisRadius, axisPt}); + histos.add("Xi/hGenXi", "hGenXi", kTH2D, {axisRadius, axisPt}); + histos.add("AntiXi/hGenAntiXi", "hGenAntiXi", kTH2D, {axisRadius, axisPt}); + histos.add("Omega/hGenOmega", "hGenOmega", kTH2D, {axisRadius, axisPt}); + histos.add("AntiOmega/hGenAntiOmega", "hGenAntiOmega", kTH2D, {axisRadius, axisPt}); } bool canDecay(const int pdgCode) @@ -103,7 +119,6 @@ struct OnTheFlyDecayer { o2::track::TrackParCov o2track; o2::upgrade::convertMCParticleToO2Track(particle, o2track, pdgDB); std::vector decayDaughters = decayer.decayParticle(pdgDB, o2track, particle.pdgCode()); - std::cout << particle.pdgCode() << " decayed into: "; for (const auto& dau : decayDaughters) { o2::track::TrackParCov dauTrack; o2::upgrade::convertOTFParticleToO2Track(dau, dauTrack, pdgDB); @@ -114,9 +129,38 @@ struct OnTheFlyDecayer { } } } - std::vector isFinalState; + + if (decayDaughters.empty()) { + LOG(error) << "Attempted to decay " << particle.pdgCode() << " but resulting vector of daugthers were empty"; + continue; + } + + + const float decayRadius2D = std::hypot(decayDaughters[0].vx, decayDaughters[0].vy); + std::cout << particle.pdgCode() << ": " << decayRadius2D << std::endl; + if (particle.pdgCode() == kK0Short) { + histos.fill(HIST("K0S/hGenK0S"), decayRadius2D, particle.pt()); + } else if (particle.pdgCode() == kLambda0) { + histos.fill(HIST("Lambda/hGenLambda"), decayRadius2D, particle.pt()); + } else if (particle.pdgCode() == kLambda0Bar) { + histos.fill(HIST("AntiLambda/hGenAntiLambda"), decayRadius2D, particle.pt()); + } else if (particle.pdgCode() == kXiMinus) { + histos.fill(HIST("Xi/hGenXi"), decayRadius2D, particle.pt()); + } else if (particle.pdgCode() == kXiPlusBar) { + histos.fill(HIST("AntiXi/hGenAntiXi"), decayRadius2D, particle.pt()); + } else if (particle.pdgCode() == kOmegaMinus) { + histos.fill(HIST("Omega/hGenOmega"), decayRadius2D, particle.pt()); + } else if (particle.pdgCode() == kOmegaPlusBar) { + histos.fill(HIST("AntiOmega/hGenAntiOmega"), decayRadius2D, particle.pt()); + } + for (size_t i = 0; i < decayDaughters.size(); ++i) { - isFinalState.push_back(!canDecay(decayDaughters[i].pdgCode)); + const o2::upgrade::OTFParticle& dau = decayDaughters[i]; + const bool isAlive = !canDecay(dau.pdgCode); + tableMcParticles(particle.globalIndex(), + isAlive, dau.pdgCode, + dau.vx, dau.vy, dau.vz, + dau.px, dau.py, dau.pz); } } } From 2d280d9cfd1524a176986a0fa784ee8db114ba6d Mon Sep 17 00:00:00 2001 From: jesgum Date: Thu, 18 Dec 2025 21:37:49 +0100 Subject: [PATCH 3/7] Add decayert to cmakelists --- ALICE3/TableProducer/OTF/CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/ALICE3/TableProducer/OTF/CMakeLists.txt b/ALICE3/TableProducer/OTF/CMakeLists.txt index b75101f977f..6fe3f41a02f 100644 --- a/ALICE3/TableProducer/OTF/CMakeLists.txt +++ b/ALICE3/TableProducer/OTF/CMakeLists.txt @@ -9,6 +9,11 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. +o2physics_add_dpl_workflow(onthefly-decayer + SOURCES onTheFlyDecayer.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2::DCAFitter O2Physics::ALICE3Core + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(onthefly-tracker SOURCES onTheFlyTracker.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2::DCAFitter O2Physics::ALICE3Core O2Physics::FastTracker From a5692451b6eaf4f63a4837e9b17e73ba86450f2f Mon Sep 17 00:00:00 2001 From: jesgum Date: Tue, 30 Dec 2025 15:11:13 +0100 Subject: [PATCH 4/7] add qa task and copy of mcparticle datamodel --- ALICE3/Core/Decayer.h | 7 +- ALICE3/Core/TrackUtilities.h | 34 ++-- ALICE3/DataModel/OTFMCParticle.h | 71 ++++++++ ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx | 174 +++++++++++++------ ALICE3/Tasks/CMakeLists.txt | 5 + ALICE3/Tasks/alice3DecayerQA.cxx | 122 +++++++++++++ 6 files changed, 342 insertions(+), 71 deletions(-) create mode 100644 ALICE3/DataModel/OTFMCParticle.h create mode 100644 ALICE3/Tasks/alice3DecayerQA.cxx diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h index efe6ecc768f..afd4e5c99ff 100644 --- a/ALICE3/Core/Decayer.h +++ b/ALICE3/Core/Decayer.h @@ -66,7 +66,7 @@ class Decayer { if (!charge) { vx = pos[0] + rxyz * (mom[0] / track.getP()); vy = pos[1] + rxyz * (mom[1] / track.getP()); - vy = pos[2] + rxyz * (mom[2] / track.getP()); + vz = pos[2] + rxyz * (mom[2] / track.getP()); px = mom[0]; py = mom[2]; } else { @@ -126,7 +126,6 @@ class Decayer { } - // Setters void setSeed(const int seed) { mRand3.SetSeed(seed); // For decay length sampling @@ -135,10 +134,6 @@ class Decayer { void setBField(const double b) { mBz = b; } - // Getters - - - private: TRandom3 mRand3; double mBz; diff --git a/ALICE3/Core/TrackUtilities.h b/ALICE3/Core/TrackUtilities.h index 9543a71bc9a..089a3553a06 100644 --- a/ALICE3/Core/TrackUtilities.h +++ b/ALICE3/Core/TrackUtilities.h @@ -39,14 +39,22 @@ void convertTLorentzVectorToO2Track(const int charge, o2::track::TrackParCov& o2track); struct OTFParticle { - int pdgCode; - float e; - float vx, vy, vz; - float px, py, pz; - - void setPDG(int pdg) { pdgCode = pdg; } - void setVxVyVz(float _vx, float _vy, float _vz) { vx = _vx; vy = _vy; vz = _vz; } - void setPxPyPzE(float _px, float _py, float _pz, float _e) { px = _px; py = _py; pz = _pz; e = _e; } + int mPdgCode; + float mE; + float mVx, mVy, mVz; + float mPx, mPy, mPz; + + void setPDG(int pdg) { mPdgCode = pdg; } + void setVxVyVz(float _vx, float _vy, float _vz) { mVx = _vx; mVy = _vy; mVz = _vz; } + void setPxPyPzE(float _px, float _py, float _pz, float _e) { mPx = _px; mPy = _py; mPz = _pz; mE = _e; } + int pdgCode() const { return mPdgCode; } + float vx() const { return mVx; } + float vy() const { return mVy; } + float vz() const { return mVz; } + float px() const { return mPx; } + float py() const { return mPy; } + float pz() const { return mPz; } + float e() const { return mE; } }; /// Function to convert a TLorentzVector into a perfect Track @@ -77,15 +85,9 @@ void convertTLorentzVectorToO2Track(int pdgCode, template void convertOTFParticleToO2Track(const OTFParticle& particle, o2::track::TrackParCov& o2track, const PdgService& pdg) { - const auto pdgInfo = pdg->GetParticle(particle.pdgCode); - int charge = 0; - if (pdgInfo != nullptr) { - charge = pdgInfo->Charge() / 3; - } - static TLorentzVector tlv; - tlv.SetPxPyPzE(particle.px, particle.py, particle.pz, particle.e); - convertTLorentzVectorToO2Track(particle.pdgCode, tlv, {particle.vx, particle.vy, particle.vz}, o2track, pdg); + tlv.SetPxPyPzE(particle.px(), particle.py(), particle.pz(), particle.e()); + convertTLorentzVectorToO2Track(particle.pdgCode(), tlv, {particle.vx(), particle.vy(), particle.vz()}, o2track, pdg); } /// Function to convert a McParticle into a perfect Track diff --git a/ALICE3/DataModel/OTFMCParticle.h b/ALICE3/DataModel/OTFMCParticle.h new file mode 100644 index 00000000000..39f932afd7d --- /dev/null +++ b/ALICE3/DataModel/OTFMCParticle.h @@ -0,0 +1,71 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file OTFMCParticle.h +/// \author Jesper Karlsson Gumprecht +/// \since 16/12/2025 +/// \brief +/// + +#ifndef ALICE3_DATAMODEL_OTFMCPARTICLE_H_ +#define ALICE3_DATAMODEL_OTFMCPARTICLE_H_ + +// O2 includes +#include "Framework/AnalysisDataModel.h" + +namespace o2::aod +{ + +namespace otfmcparticle +{ + DECLARE_SOA_COLUMN(Phi, phi, float); + DECLARE_SOA_COLUMN(Eta, eta, float); + DECLARE_SOA_COLUMN(Pt, pt, float); + DECLARE_SOA_COLUMN(P, p, float); + DECLARE_SOA_COLUMN(Y, y, float); +} + +DECLARE_SOA_TABLE_FULL(McParticlesWithDau, "McParticlesWithDau", "AOD", "MCPARTICLEWITHDAU", + o2::soa::Index<>, + mcparticle::McCollisionId, + mcparticle::PdgCode, + mcparticle::StatusCode, + mcparticle::Flags, + mcparticle::MothersIds, + mcparticle::DaughtersIdSlice, + mcparticle::Weight, + mcparticle::Px, + mcparticle::Py, + mcparticle::Pz, + mcparticle::E, + mcparticle::Vx, + mcparticle::Vy, + mcparticle::Vz, + mcparticle::Vt, + otfmcparticle::Phi, + otfmcparticle::Eta, + otfmcparticle::Pt, + otfmcparticle::P, + otfmcparticle::Y, + mcparticle::PVector, + mcparticle::ProducedByGenerator, + mcparticle::FromBackgroundEvent, + mcparticle::GetGenStatusCode, + mcparticle::GetHepMCStatusCode, + mcparticle::GetProcess, + mcparticle::IsPhysicalPrimary); + +using McParticleWithDau = McParticlesWithDau::iterator; + +} // namespace o2::aod + +#endif // ALICE3_DATAMODEL_OTFMCPARTICLE_H_ \ No newline at end of file diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx index 391e2096a82..18ef0bd86d5 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -40,6 +40,8 @@ #include #include +#include +#include #include using namespace o2; @@ -65,11 +67,11 @@ static const std::vector pdgCodes{kK0Short, kOmegaPlusBar}; struct OnTheFlyDecayer { - Produces tableMcParticles; + Produces tableMcParticlesWithDau; o2::upgrade::Decayer decayer; - Service ccdb; Service pdgDB; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + std::map> mDecayDaughters; Configurable seed{"seed", 0, "Set seed for particle decayer"}; Configurable magneticField{"magneticField", 20., "Magnetic field (kG)"}; @@ -78,15 +80,12 @@ struct OnTheFlyDecayer { "Enable option for particle to be decayed: 0 - no, 1 - yes"}; - ConfigurableAxis axisRadius{"axisRadius", {1000, 0, 100}, "Radius"}; 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"}; std::vector mEnabledDecays; void init(o2::framework::InitContext&) { - ccdb->setURL("http://alice-ccdb.cern.ch"); - ccdb->setTimestamp(-1); decayer.setSeed(seed); decayer.setBField(magneticField); for (int i = 0; i < kNumDecays; ++i) { @@ -95,13 +94,13 @@ struct OnTheFlyDecayer { } } - histos.add("K0S/hGenK0S", "hGenK0S", kTH2D, {axisRadius, axisPt}); - histos.add("Lambda/hGenLambda", "hGenLambda", kTH2D, {axisRadius, axisPt}); - histos.add("AntiLambda/hGenAntiLambda", "hGenAntiLambda", kTH2D, {axisRadius, axisPt}); - histos.add("Xi/hGenXi", "hGenXi", kTH2D, {axisRadius, axisPt}); - histos.add("AntiXi/hGenAntiXi", "hGenAntiXi", kTH2D, {axisRadius, axisPt}); - histos.add("Omega/hGenOmega", "hGenOmega", kTH2D, {axisRadius, axisPt}); - histos.add("AntiOmega/hGenAntiOmega", "hGenAntiOmega", kTH2D, {axisRadius, axisPt}); + histos.add("K0S/hGenK0S", "hGenK0S", kTH1D, {axisPt}); + histos.add("Lambda/hGenLambda", "hGenLambda", kTH1D, {axisPt}); + histos.add("AntiLambda/hGenAntiLambda", "hGenAntiLambda", kTH1D, {axisPt}); + histos.add("Xi/hGenXi", "hGenXi", kTH1D, {axisPt}); + histos.add("AntiXi/hGenAntiXi", "hGenAntiXi", kTH1D, {axisPt}); + histos.add("Omega/hGenOmega", "hGenOmega", kTH1D, {axisPt}); + histos.add("AntiOmega/hGenAntiOmega", "hGenAntiOmega", kTH1D, {axisPt}); } bool canDecay(const int pdgCode) @@ -111,56 +110,133 @@ struct OnTheFlyDecayer { void process(aod::McCollision const&, aod::McParticles const& mcParticles) { - for (const auto& particle : mcParticles) { - if (!canDecay(particle.pdgCode())) { - continue; - } - - o2::track::TrackParCov o2track; - o2::upgrade::convertMCParticleToO2Track(particle, o2track, pdgDB); - std::vector decayDaughters = decayer.decayParticle(pdgDB, o2track, particle.pdgCode()); - for (const auto& dau : decayDaughters) { - o2::track::TrackParCov dauTrack; - o2::upgrade::convertOTFParticleToO2Track(dau, dauTrack, pdgDB); - if (canDecay(dau.pdgCode)) { - std::vector cascadingDaughers = decayer.decayParticle(pdgDB, dauTrack, dau.pdgCode); - for (const auto& daudau : cascadingDaughers) { - decayDaughters.push_back(daudau); + mDecayDaughters.clear(); + u_int64_t nStoredDaughters = 0; + for (int64_t index{0}; index < mcParticles.size(); ++index) { + const auto& particle = mcParticles.iteratorAt(index); + std::vector decayDaughters; + static constexpr int kMaxNestedDecays = 10; + int nDecays = 0; + if (canDecay(particle.pdgCode())) { + o2::track::TrackParCov o2track; + o2::upgrade::convertMCParticleToO2Track(particle, o2track, pdgDB); + decayDaughters = decayer.decayParticle(pdgDB, o2track, particle.pdgCode()); + for (size_t idau{0}; idau < decayDaughters.size(); ++idau) { + const auto& dau = decayDaughters[idau]; + o2::track::TrackParCov dauTrack; + o2::upgrade::convertOTFParticleToO2Track(dau, dauTrack, pdgDB); + if (canDecay(dau.pdgCode())) { + std::vector cascadingDaughers = decayer.decayParticle(pdgDB, dauTrack, dau.pdgCode()); + for (const auto& daudau : cascadingDaughers) { + decayDaughters.push_back(daudau); + if (kMaxNestedDecays < ++nDecays) { + LOG(error) << "Seemingly stuck trying to perpetually decay products from pdg: " << particle.pdgCode(); + } + } } } - } - if (decayDaughters.empty()) { - LOG(error) << "Attempted to decay " << particle.pdgCode() << " but resulting vector of daugthers were empty"; - continue; + if (decayDaughters.empty()) { + LOG(error) << "Attempted to decay " << particle.pdgCode() << " but resulting vector of daugthers were empty"; + continue; + } } - const float decayRadius2D = std::hypot(decayDaughters[0].vx, decayDaughters[0].vy); - std::cout << particle.pdgCode() << ": " << decayRadius2D << std::endl; if (particle.pdgCode() == kK0Short) { - histos.fill(HIST("K0S/hGenK0S"), decayRadius2D, particle.pt()); + histos.fill(HIST("K0S/hGenK0S"), particle.pt()); } else if (particle.pdgCode() == kLambda0) { - histos.fill(HIST("Lambda/hGenLambda"), decayRadius2D, particle.pt()); + histos.fill(HIST("Lambda/hGenLambda"), particle.pt()); } else if (particle.pdgCode() == kLambda0Bar) { - histos.fill(HIST("AntiLambda/hGenAntiLambda"), decayRadius2D, particle.pt()); + histos.fill(HIST("AntiLambda/hGenAntiLambda"), particle.pt()); } else if (particle.pdgCode() == kXiMinus) { - histos.fill(HIST("Xi/hGenXi"), decayRadius2D, particle.pt()); + histos.fill(HIST("Xi/hGenXi"), particle.pt()); } else if (particle.pdgCode() == kXiPlusBar) { - histos.fill(HIST("AntiXi/hGenAntiXi"), decayRadius2D, particle.pt()); + histos.fill(HIST("AntiXi/hGenAntiXi"), particle.pt()); } else if (particle.pdgCode() == kOmegaMinus) { - histos.fill(HIST("Omega/hGenOmega"), decayRadius2D, particle.pt()); + histos.fill(HIST("Omega/hGenOmega"), particle.pt()); } else if (particle.pdgCode() == kOmegaPlusBar) { - histos.fill(HIST("AntiOmega/hGenAntiOmega"), decayRadius2D, particle.pt()); + histos.fill(HIST("AntiOmega/hGenAntiOmega"), particle.pt()); } - - for (size_t i = 0; i < decayDaughters.size(); ++i) { - const o2::upgrade::OTFParticle& dau = decayDaughters[i]; - const bool isAlive = !canDecay(dau.pdgCode); - tableMcParticles(particle.globalIndex(), - isAlive, dau.pdgCode, - dau.vx, dau.vy, dau.vz, - dau.px, dau.py, dau.pz); + + int daughtersIdSlice[2]; + if (canDecay(particle.pdgCode())) { + daughtersIdSlice[0] = static_cast(mcParticles.size() + nStoredDaughters); + daughtersIdSlice[1] = static_cast(mcParticles.size() + nStoredDaughters + decayDaughters.size()); + } else { + daughtersIdSlice[0] = static_cast(particle.daughtersIds()[0]); + 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()); + float eta; // Conditional 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()); + float y; // Conditional as https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943 + + if ((p - particle.pz()) < 1e-7f) { + eta = (particle.pz() < 0.0f) ? -100.0f : 100.0f; + } else { + eta = 0.5f * std::log((p + particle.pz()) / (p - particle.pz())); + } + + if ((particle.e() - particle.pz()) < 1e-7f) { + y = (particle.pz() < 0.0f) ? -100.0f : 100.0f; + } else { + y = 0.5f * std::log((particle.e() + particle.pz()) / (particle.e() - particle.pz())); + } + + // 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); + } + + 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(); + 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()); + 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()); + float y; // Conditional as https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943 + + if ((p - dau.pz()) < 1e-7f) { + eta = (dau.pz() < 0.0f) ? -100.0f : 100.0f; + } else { + eta = 0.5f * std::log((p + dau.pz()) / (p - dau.pz())); + } + + if ((dau.e() - dau.pz()) < 1e-7f) { + y = (dau.pz() < 0.0f) ? -100.0f : 100.0f; + } else { + y = 0.5f * std::log((dau.e() + dau.pz()) / (dau.e() - dau.pz())); + } + + + // TODO: Particle status code + // TODO: Expression columns + 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); } } } @@ -170,4 +246,4 @@ struct OnTheFlyDecayer { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; -} +} \ No newline at end of file diff --git a/ALICE3/Tasks/CMakeLists.txt b/ALICE3/Tasks/CMakeLists.txt index 56c162a1f97..f1158275ec8 100644 --- a/ALICE3/Tasks/CMakeLists.txt +++ b/ALICE3/Tasks/CMakeLists.txt @@ -88,3 +88,8 @@ o2physics_add_dpl_workflow(alice3-strangeness SOURCES alice3-strangeness.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(alice3-decayer-qa + SOURCES alice3DecayerQA.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/ALICE3/Tasks/alice3DecayerQA.cxx b/ALICE3/Tasks/alice3DecayerQA.cxx new file mode 100644 index 00000000000..04c133c3e05 --- /dev/null +++ b/ALICE3/Tasks/alice3DecayerQA.cxx @@ -0,0 +1,122 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file alice3DecayerQA.cxx +/// +/// \brief QA task for otf decayer +/// +/// \author Jesper Karlsson Gumprecht +/// \since Dec 23, 2025 +/// + +#include +#include + +#include +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/ASoAHelpers.h" + +#include "ALICE3/DataModel/OTFMCParticle.h" + +#include + + +using namespace o2; +using namespace o2::framework; + +struct Alice3DecayerQA { + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.101f, 1.131f}, ""}; + ConfigurableAxis axisXiMass{"axisXiMass", {200, 1.22f, 1.42f}, ""}; + ConfigurableAxis axisRadius{"axisRadius", {1000, 0, 100}, "Radius"}; + 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"}; + + Partition trueEl = aod::mcparticle::pdgCode == static_cast(kElectron); + Partition trueMu = aod::mcparticle::pdgCode == static_cast(kMuonMinus); + Partition truePi = aod::mcparticle::pdgCode == static_cast(kPiPlus); + Partition trueKa = aod::mcparticle::pdgCode == static_cast(kKMinus); + Partition truePr = aod::mcparticle::pdgCode == static_cast(kProton); + + Partition trueElWithDau = aod::mcparticle::pdgCode == static_cast(kElectron); + Partition trueMuWithDau = aod::mcparticle::pdgCode == static_cast(kMuonMinus); + Partition truePiWithDau = aod::mcparticle::pdgCode == static_cast(kPiPlus); + Partition trueKaWithDau = aod::mcparticle::pdgCode == static_cast(kKMinus); + Partition truePrWithDau = aod::mcparticle::pdgCode == static_cast(kProton); + + + 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("Lambda/hGenLambda", "hGenLambda", kTH2D, {axisRadius, axisPt}); + histos.add("Xi/hGenXi", "hGenXi", kTH2D, {axisRadius, axisPt}); + } + + void processMC(const aod::McParticles&) + { + for (auto const& particle : trueEl) { + histos.fill(HIST("DefaultMC/hElPt"), particle.pt()); + } + for (auto const& particle : trueMu) { + histos.fill(HIST("DefaultMC/hMuPt"), particle.pt()); + } + for (auto const& particle : truePi) { + histos.fill(HIST("DefaultMC/hPiPt"), particle.pt()); + } + for (auto const& particle : trueKa) { + histos.fill(HIST("DefaultMC/hKaPt"), particle.pt()); + } + for (auto const& particle : truePr) { + histos.fill(HIST("DefaultMC/hPrPt"), particle.pt()); + } + } + + void processMCWithDau(const aod::McParticlesWithDau&) + { + for (auto const& particle : trueElWithDau) { + histos.fill(HIST("MCWithDau/hElPt"), particle.pt()); + } + for (auto const& particle : trueMuWithDau) { + histos.fill(HIST("MCWithDau/hMuPt"), particle.pt()); + } + for (auto const& particle : truePiWithDau) { + histos.fill(HIST("MCWithDau/hPiPt"), particle.pt()); + } + for (auto const& particle : trueKaWithDau) { + histos.fill(HIST("MCWithDau/hKaPt"), particle.pt()); + } + for (auto const& particle : truePrWithDau) { + histos.fill(HIST("MCWithDau/hPrPt"), particle.pt()); + } + } + + PROCESS_SWITCH(Alice3DecayerQA, processMC, "fill MC-only histograms", true); + PROCESS_SWITCH(Alice3DecayerQA, processMCWithDau, "fill MC-only histograms", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& ctx) +{ + return WorkflowSpec{adaptAnalysisTask(ctx)}; +} From 2768d7d40e308057cdd8c85dc5a48691e07ca65d Mon Sep 17 00:00:00 2001 From: jesgum Date: Tue, 30 Dec 2025 15:27:45 +0100 Subject: [PATCH 5/7] Formatting --- ALICE3/Core/Decayer.h | 22 +++++----- ALICE3/Core/TrackUtilities.h | 42 +++++++++++++------- ALICE3/DataModel/OTFMCParticle.h | 38 +++++++++--------- ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx | 31 +++++++-------- ALICE3/Tasks/alice3DecayerQA.cxx | 29 +++++--------- 5 files changed, 82 insertions(+), 80 deletions(-) diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h index afd4e5c99ff..cb5fdf0cb26 100644 --- a/ALICE3/Core/Decayer.h +++ b/ALICE3/Core/Decayer.h @@ -19,6 +19,7 @@ #define ALICE3_CORE_DECAYER_H_ #include "ALICE3/Core/TrackUtilities.h" + #include "ReconstructionDataFormats/Track.h" #include @@ -28,19 +29,19 @@ #include #include +#include +#include #include #include #include -#include -#include - namespace o2 { namespace upgrade { -class Decayer { +class Decayer +{ public: // Default constructor Decayer() = default; @@ -75,11 +76,11 @@ class Decayer { track.getCircleParams(mBz, circle, sna, csa); const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl()); const double theta = rxy / circle.rC; - + vx = ((pos[0] - circle.xC) * std::cos(theta) - (pos[1] - circle.yC) * std::sin(theta)) + circle.xC; vy = ((pos[1] - circle.yC) * std::cos(theta) + (pos[0] - circle.xC) * std::sin(theta)) + circle.yC; vz = mom[2] + rxyz * (mom[2] / track.getP()); - + px = mom[0] * std::cos(theta) - mom[1] * std::sin(theta); py = mom[1] * std::cos(theta) + mom[0] * std::sin(theta); } @@ -89,7 +90,7 @@ class Decayer { for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) { brTotal += particleInfo->DecayChannel(ch)->BranchingRatio(); } - + double brSum = 0.; std::vector dauMasses; std::vector pdgCodesDaughters; @@ -125,21 +126,20 @@ class Decayer { return decayProducts; } - void setSeed(const int seed) { - mRand3.SetSeed(seed); // For decay length sampling + mRand3.SetSeed(seed); // For decay length sampling gRandom->SetSeed(seed); // For TGenPhaseSpace } - void setBField(const double b) { mBz = b; } + void setBField(const double b) { mBz = b; } private: TRandom3 mRand3; double mBz; }; -} // namespace fastsim +} // namespace upgrade } // namespace o2 #endif // ALICE3_CORE_DECAYER_H_ diff --git a/ALICE3/Core/TrackUtilities.h b/ALICE3/Core/TrackUtilities.h index 089a3553a06..dcfb3bfb245 100644 --- a/ALICE3/Core/TrackUtilities.h +++ b/ALICE3/Core/TrackUtilities.h @@ -21,32 +21,36 @@ #include "ReconstructionDataFormats/Track.h" #include "TLorentzVector.h" -#include - +#include namespace o2::upgrade { -/// Function to convert a TLorentzVector into a perfect Track -/// \param charge particle charge (integer) -/// \param particle the particle to convert (TLorentzVector) -/// \param productionVertex where the particle was produced -/// \param o2track the address of the resulting TrackParCov -void convertTLorentzVectorToO2Track(const int charge, - const TLorentzVector particle, - const std::vector productionVertex, - o2::track::TrackParCov& o2track); - +/// Struct to store mc info for the otf decayer struct OTFParticle { int mPdgCode; float mE; float mVx, mVy, mVz; float mPx, mPy, mPz; + // Setters void setPDG(int pdg) { mPdgCode = pdg; } - void setVxVyVz(float _vx, float _vy, float _vz) { mVx = _vx; mVy = _vy; mVz = _vz; } - void setPxPyPzE(float _px, float _py, float _pz, float _e) { mPx = _px; mPy = _py; mPz = _pz; mE = _e; } + void setVxVyVz(float vx, float vy, float vz) + { + mVx = vx; + mVy = vy; + mVz = vz; + } + void setPxPyPzE(float px, float py, float pz, float e) + { + mPx = px; + mPy = py; + mPz = pz; + mE = e; + } + + // Getters int pdgCode() const { return mPdgCode; } float vx() const { return mVx; } float vy() const { return mVy; } @@ -57,6 +61,16 @@ struct OTFParticle { float e() const { return mE; } }; +/// Function to convert a TLorentzVector into a perfect Track +/// \param charge particle charge (integer) +/// \param particle the particle to convert (TLorentzVector) +/// \param productionVertex where the particle was produced +/// \param o2track the address of the resulting TrackParCov +void convertTLorentzVectorToO2Track(const int charge, + const TLorentzVector particle, + const std::vector productionVertex, + o2::track::TrackParCov& o2track); + /// Function to convert a TLorentzVector into a perfect Track /// \param pdgCode particle pdg /// \param particle the particle to convert (TLorentzVector) diff --git a/ALICE3/DataModel/OTFMCParticle.h b/ALICE3/DataModel/OTFMCParticle.h index 39f932afd7d..b95a394473e 100644 --- a/ALICE3/DataModel/OTFMCParticle.h +++ b/ALICE3/DataModel/OTFMCParticle.h @@ -13,7 +13,7 @@ /// \file OTFMCParticle.h /// \author Jesper Karlsson Gumprecht /// \since 16/12/2025 -/// \brief +/// \brief /// #ifndef ALICE3_DATAMODEL_OTFMCPARTICLE_H_ @@ -27,29 +27,29 @@ namespace o2::aod namespace otfmcparticle { - DECLARE_SOA_COLUMN(Phi, phi, float); - DECLARE_SOA_COLUMN(Eta, eta, float); - DECLARE_SOA_COLUMN(Pt, pt, float); - DECLARE_SOA_COLUMN(P, p, float); - DECLARE_SOA_COLUMN(Y, y, float); -} +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(P, p, float); +DECLARE_SOA_COLUMN(Y, y, float); +} // namespace otfmcparticle DECLARE_SOA_TABLE_FULL(McParticlesWithDau, "McParticlesWithDau", "AOD", "MCPARTICLEWITHDAU", - o2::soa::Index<>, + o2::soa::Index<>, mcparticle::McCollisionId, - mcparticle::PdgCode, - mcparticle::StatusCode, + mcparticle::PdgCode, + mcparticle::StatusCode, mcparticle::Flags, - mcparticle::MothersIds, - mcparticle::DaughtersIdSlice, + mcparticle::MothersIds, + mcparticle::DaughtersIdSlice, mcparticle::Weight, - mcparticle::Px, - mcparticle::Py, - mcparticle::Pz, + mcparticle::Px, + mcparticle::Py, + mcparticle::Pz, mcparticle::E, - mcparticle::Vx, - mcparticle::Vy, - mcparticle::Vz, + mcparticle::Vx, + mcparticle::Vy, + mcparticle::Vz, mcparticle::Vt, otfmcparticle::Phi, otfmcparticle::Eta, @@ -68,4 +68,4 @@ using McParticleWithDau = McParticlesWithDau::iterator; } // namespace o2::aod -#endif // ALICE3_DATAMODEL_OTFMCPARTICLE_H_ \ No newline at end of file +#endif // ALICE3_DATAMODEL_OTFMCPARTICLE_H_ diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx index 18ef0bd86d5..2d8af38fa06 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -39,9 +39,9 @@ #include -#include #include #include +#include #include using namespace o2; @@ -68,25 +68,24 @@ static const std::vector pdgCodes{kK0Short, struct OnTheFlyDecayer { Produces tableMcParticlesWithDau; + o2::upgrade::Decayer decayer; Service pdgDB; - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; std::map> mDecayDaughters; Configurable seed{"seed", 0, "Set seed for particle decayer"}; Configurable magneticField{"magneticField", 20., "Magnetic field (kG)"}; Configurable> enabledDecays{"enabledDecays", - {defaultParameters[0], kNumDecays, kNumParameters, particleNames, parameterNames}, - "Enable option for particle to be decayed: 0 - no, 1 - yes"}; - + {defaultParameters[0], kNumDecays, kNumParameters, particleNames, parameterNames}, + "Enable option for particle to be decayed: 0 - no, 1 - yes"}; + 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"}; - std::vector mEnabledDecays; void init(o2::framework::InitContext&) { - decayer.setSeed(seed); + decayer.setSeed(seed); decayer.setBField(magneticField); for (int i = 0; i < kNumDecays; ++i) { if (enabledDecays->get(particleNames[i].c_str(), "enable")) { @@ -106,7 +105,7 @@ struct OnTheFlyDecayer { bool canDecay(const int pdgCode) { return std::find(mEnabledDecays.begin(), mEnabledDecays.end(), pdgCode) != mEnabledDecays.end(); - } + } void process(aod::McCollision const&, aod::McParticles const& mcParticles) { @@ -141,8 +140,7 @@ struct OnTheFlyDecayer { continue; } } - - + if (particle.pdgCode() == kK0Short) { histos.fill(HIST("K0S/hGenK0S"), particle.pt()); } else if (particle.pdgCode() == kLambda0) { @@ -192,9 +190,9 @@ struct OnTheFlyDecayer { // TODO: Particle status code // TODO: Expression columns - tableMcParticlesWithDau(particle.mcCollisionId(), particle.pdgCode(), particle.statusCode(), + tableMcParticlesWithDau(particle.mcCollisionId(), particle.pdgCode(), particle.statusCode(), particle.flags(), motherSpan, daughtersIdSlice, particle.weight(), - particle.px(), particle.py(), particle.pz(), particle.e(), + particle.px(), particle.py(), particle.pz(), particle.e(), particle.vx(), particle.vy(), particle.vz(), particle.vt(), phi, eta, pt, p, y); } @@ -208,7 +206,7 @@ struct OnTheFlyDecayer { } auto mother = mcParticles.iteratorAt(index); - std::vector motherIds = { static_cast(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()); @@ -229,11 +227,11 @@ struct OnTheFlyDecayer { y = 0.5f * std::log((dau.e() + dau.pz()) / (dau.e() - dau.pz())); } - // TODO: Particle status code // TODO: Expression columns + // TODO: vt tableMcParticlesWithDau(mother.mcCollisionId(), dau.pdgCode(), 1, - mother.flags(), motherSpan, daughtersIdSlice, mother.weight(), + 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); @@ -242,8 +240,7 @@ struct OnTheFlyDecayer { } }; - WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; -} \ No newline at end of file +} diff --git a/ALICE3/Tasks/alice3DecayerQA.cxx b/ALICE3/Tasks/alice3DecayerQA.cxx index 04c133c3e05..7eff59fe2d0 100644 --- a/ALICE3/Tasks/alice3DecayerQA.cxx +++ b/ALICE3/Tasks/alice3DecayerQA.cxx @@ -13,34 +13,29 @@ /// /// \brief QA task for otf decayer /// -/// \author Jesper Karlsson Gumprecht +/// \author Jesper Karlsson Gumprecht /// \since Dec 23, 2025 /// -#include -#include +#include "ALICE3/DataModel/OTFMCParticle.h" -#include +#include "Framework/ASoAHelpers.h" #include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "Framework/HistogramRegistry.h" #include "Framework/ConfigParamRegistry.h" -#include "Framework/ASoAHelpers.h" - -#include "ALICE3/DataModel/OTFMCParticle.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" +#include #include +#include +#include using namespace o2; using namespace o2::framework; struct Alice3DecayerQA { HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - - ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.101f, 1.131f}, ""}; - ConfigurableAxis axisXiMass{"axisXiMass", {200, 1.22f, 1.42f}, ""}; - ConfigurableAxis axisRadius{"axisRadius", {1000, 0, 100}, "Radius"}; 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"}; Partition trueEl = aod::mcparticle::pdgCode == static_cast(kElectron); @@ -48,14 +43,13 @@ struct Alice3DecayerQA { Partition truePi = aod::mcparticle::pdgCode == static_cast(kPiPlus); Partition trueKa = aod::mcparticle::pdgCode == static_cast(kKMinus); Partition truePr = aod::mcparticle::pdgCode == static_cast(kProton); - + Partition trueElWithDau = aod::mcparticle::pdgCode == static_cast(kElectron); Partition trueMuWithDau = aod::mcparticle::pdgCode == static_cast(kMuonMinus); Partition truePiWithDau = aod::mcparticle::pdgCode == static_cast(kPiPlus); Partition trueKaWithDau = aod::mcparticle::pdgCode == static_cast(kKMinus); Partition truePrWithDau = aod::mcparticle::pdgCode == static_cast(kProton); - void init(o2::framework::InitContext&) { histos.add("DefaultMC/hElPt", "hElPt", kTH1D, {axisPt}); @@ -69,9 +63,6 @@ struct Alice3DecayerQA { histos.add("MCWithDau/hPiPt", "hPiPt", kTH1D, {axisPt}); histos.add("MCWithDau/hKaPt", "hKaPt", kTH1D, {axisPt}); histos.add("MCWithDau/hPrPt", "hPrPt", kTH1D, {axisPt}); - - histos.add("Lambda/hGenLambda", "hGenLambda", kTH2D, {axisRadius, axisPt}); - histos.add("Xi/hGenXi", "hGenXi", kTH2D, {axisRadius, axisPt}); } void processMC(const aod::McParticles&) @@ -113,7 +104,7 @@ struct Alice3DecayerQA { } PROCESS_SWITCH(Alice3DecayerQA, processMC, "fill MC-only histograms", true); - PROCESS_SWITCH(Alice3DecayerQA, processMCWithDau, "fill MC-only histograms", true); + PROCESS_SWITCH(Alice3DecayerQA, processMCWithDau, "fill MC-with-dau histograms", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& ctx) From e45f1a80c8072c1764851e9957e93a9c396b5fa8 Mon Sep 17 00:00:00 2001 From: jesgum Date: Wed, 7 Jan 2026 16:06:46 +0100 Subject: [PATCH 6/7] cleanup --- ALICE3/Core/Decayer.h | 1 + ALICE3/DataModel/OTFMCParticle.h | 2 +- ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx | 27 ++++++++++---------- 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h index cb5fdf0cb26..c6358aa5c1c 100644 --- a/ALICE3/Core/Decayer.h +++ b/ALICE3/Core/Decayer.h @@ -13,6 +13,7 @@ /// \file Decayer.h /// \author Jesper Karlsson Gumprecht /// \since 15/12/2025 +/// \brief Basic class to handle short-lived particle decays in the fast simulation /// #ifndef ALICE3_CORE_DECAYER_H_ diff --git a/ALICE3/DataModel/OTFMCParticle.h b/ALICE3/DataModel/OTFMCParticle.h index b95a394473e..92f5af37c3e 100644 --- a/ALICE3/DataModel/OTFMCParticle.h +++ b/ALICE3/DataModel/OTFMCParticle.h @@ -13,7 +13,7 @@ /// \file OTFMCParticle.h /// \author Jesper Karlsson Gumprecht /// \since 16/12/2025 -/// \brief +/// \brief Redefinition of the mcparticles table specifically for the fast sim /// #ifndef ALICE3_DATAMODEL_OTFMCPARTICLE_H_ diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx index 2d8af38fa06..8f1e1448d0a 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -47,9 +47,10 @@ using namespace o2; using namespace o2::framework; -static constexpr int kNumDecays = 7; -static constexpr int kNumParameters = 1; -static constexpr int defaultParameters[kNumDecays][kNumParameters]{{1}, {1}, {1}, {1}, {1}, {1}, {1}}; +static constexpr int NumDecays = 7; +static constexpr int NumParameters = 1; +static constexpr int DefaultParameters[NumDecays][NumParameters]{{1}, {1}, {1}, {1}, {1}, {1}, {1}}; +static constexpr float VerySmall = 1e-7f; static const std::vector parameterNames{"enable"}; static const std::vector particleNames{"K0s", "Lambda", @@ -76,7 +77,7 @@ struct OnTheFlyDecayer { Configurable seed{"seed", 0, "Set seed for particle decayer"}; Configurable magneticField{"magneticField", 20., "Magnetic field (kG)"}; Configurable> enabledDecays{"enabledDecays", - {defaultParameters[0], kNumDecays, kNumParameters, particleNames, parameterNames}, + {DefaultParameters[0], NumDecays, NumParameters, particleNames, parameterNames}, "Enable option for particle to be decayed: 0 - no, 1 - yes"}; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -87,7 +88,7 @@ struct OnTheFlyDecayer { { decayer.setSeed(seed); decayer.setBField(magneticField); - for (int i = 0; i < kNumDecays; ++i) { + for (int i = 0; i < NumDecays; ++i) { if (enabledDecays->get(particleNames[i].c_str(), "enable")) { mEnabledDecays.push_back(pdgCodes[i]); } @@ -114,7 +115,7 @@ struct OnTheFlyDecayer { for (int64_t index{0}; index < mcParticles.size(); ++index) { const auto& particle = mcParticles.iteratorAt(index); std::vector decayDaughters; - static constexpr int kMaxNestedDecays = 10; + static constexpr int MaxNestedDecays = 10; int nDecays = 0; if (canDecay(particle.pdgCode())) { o2::track::TrackParCov o2track; @@ -128,7 +129,7 @@ struct OnTheFlyDecayer { std::vector cascadingDaughers = decayer.decayParticle(pdgDB, dauTrack, dau.pdgCode()); for (const auto& daudau : cascadingDaughers) { decayDaughters.push_back(daudau); - if (kMaxNestedDecays < ++nDecays) { + if (MaxNestedDecays < ++nDecays) { LOG(error) << "Seemingly stuck trying to perpetually decay products from pdg: " << particle.pdgCode(); } } @@ -171,18 +172,18 @@ struct OnTheFlyDecayer { nStoredDaughters += decayDaughters.size(); float phi = o2::constants::math::PI + std::atan2(-1.0f * particle.py(), -1.0f * particle.px()); - float eta; // Conditional as https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1922 + 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()); - float y; // Conditional as https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943 + float y; // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943 - if ((p - particle.pz()) < 1e-7f) { + if ((p - particle.pz()) < VerySmall) { eta = (particle.pz() < 0.0f) ? -100.0f : 100.0f; } else { eta = 0.5f * std::log((p + particle.pz()) / (p - particle.pz())); } - if ((particle.e() - particle.pz()) < 1e-7f) { + if ((particle.e() - particle.pz()) < VerySmall) { y = (particle.pz() < 0.0f) ? -100.0f : 100.0f; } else { y = 0.5f * std::log((particle.e() + particle.pz()) / (particle.e() - particle.pz())); @@ -215,13 +216,13 @@ struct OnTheFlyDecayer { 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()) < 1e-7f) { + if ((p - dau.pz()) < VerySmall) { eta = (dau.pz() < 0.0f) ? -100.0f : 100.0f; } else { eta = 0.5f * std::log((p + dau.pz()) / (p - dau.pz())); } - if ((dau.e() - dau.pz()) < 1e-7f) { + if ((dau.e() - dau.pz()) < VerySmall) { y = (dau.pz() < 0.0f) ? -100.0f : 100.0f; } else { y = 0.5f * std::log((dau.e() + dau.pz()) / (dau.e() - dau.pz())); From 786dc3382d9f8202247f1ebafa8686e53bd9ba4f Mon Sep 17 00:00:00 2001 From: jesgum Date: Fri, 9 Jan 2026 11:00:20 +0100 Subject: [PATCH 7/7] megalinter fix --- ALICE3/Core/Decayer.h | 1 - ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h index c6358aa5c1c..22cc46139cd 100644 --- a/ALICE3/Core/Decayer.h +++ b/ALICE3/Core/Decayer.h @@ -33,7 +33,6 @@ #include #include #include -#include #include namespace o2 diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx index 8f1e1448d0a..f61bab5605f 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -40,6 +40,7 @@ #include #include +#include #include #include #include