add qa task and copy of mcparticle datamodel

jesgum · jesgum · commit a5692451b6ea · 2026-01-07T15:40:43.000+01:00
diff --git 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
@@ -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 <typename PdgService>
 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
@@ -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::Px, mcparticle::Py, mcparticle::Pz>,
+                       mcparticle::ProducedByGenerator<mcparticle::Flags>,
+                       mcparticle::FromBackgroundEvent<mcparticle::Flags>,
+                       mcparticle::GetGenStatusCode<mcparticle::Flags, mcparticle::StatusCode>,
+                       mcparticle::GetHepMCStatusCode<mcparticle::Flags, mcparticle::StatusCode>,
+                       mcparticle::GetProcess<mcparticle::Flags, mcparticle::StatusCode>,
+                       mcparticle::IsPhysicalPrimary<mcparticle::Flags>);
+
+using McParticleWithDau = McParticlesWithDau::iterator;
+
+} // namespace o2::aod
+
+#endif // ALICE3_DATAMODEL_OTFMCPARTICLE_H_
diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx
@@ -40,6 +40,8 @@
 #include <TPDGCode.h>
 
 #include <string>
+#include <gsl/span>
+#include <span>
 #include <vector>
 
 using namespace o2;
@@ -65,11 +67,11 @@ static const std::vector<int> pdgCodes{kK0Short,
                                        kOmegaPlusBar};
 
 struct OnTheFlyDecayer {
-  Produces<aod::OTFMCParticles> tableMcParticles;
+  Produces<aod::McParticlesWithDau> tableMcParticlesWithDau;
   o2::upgrade::Decayer decayer;
-  Service<o2::ccdb::BasicCCDBManager> ccdb;
   Service<o2::framework::O2DatabasePDG> pdgDB;
   HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
+  std::map<int64_t, std::vector<o2::upgrade::OTFParticle>> mDecayDaughters;
 
   Configurable<int> seed{"seed", 0, "Set seed for particle decayer"};
   Configurable<float> 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<int> 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<o2::upgrade::OTFParticle> 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<o2::upgrade::OTFParticle> 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<o2::upgrade::OTFParticle> 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<o2::upgrade::OTFParticle> 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<int>(mcParticles.size() + nStoredDaughters);
+        daughtersIdSlice[1] = static_cast<int>(mcParticles.size() + nStoredDaughters + decayDaughters.size());
+      } else {
+        daughtersIdSlice[0] = static_cast<int>(particle.daughtersIds()[0]);
+        daughtersIdSlice[1] = static_cast<int>(particle.daughtersIds()[1]);
+      }
+
+      std::span<const int> 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<int> motherIds = { static_cast<int>(index) };
+        std::span<const int> 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<OnTheFlyDecayer>(cfgc)};
-}
+}
diff --git 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
