[ALICE3] Add procces function to onTheFlyTracker.cxx for mcParticle with daughters

ALICE3/Core/Decayer.h

@@ -23,10 +23,10 @@
 
 #include "ReconstructionDataFormats/Track.h"
 
 #include <TDatabasePDG.h>
 #include <TDecayChannel.h>
 #include <TGenPhaseSpace.h>
 #include <TLorentzVector.h>
 #include <TParticlePDG.h>
 #include <TRandom3.h>
 
@@ -69,7 +69,7 @@
       vy = pos[1] + rxyz * (mom[1] / track.getP());
       vz = pos[2] + rxyz * (mom[2] / track.getP());
       px = mom[0];
-      py = mom[2];
+      py = mom[1];
     } else {
       float sna, csa;
       o2::math_utils::CircleXYf_t circle;
@@ -108,7 +108,7 @@
       }
     }
 
     TLorentzVector tlv(px, py, mom[2], e);
     TGenPhaseSpace decay;
     decay.SetDecay(tlv, dauMasses.size(), dauMasses.data());
     decay.Generate();
@@ -116,7 +116,7 @@
     std::vector<o2::upgrade::OTFParticle> 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());

ALICE3/Core/TrackUtilities.h

@@ -20,7 +20,7 @@
 
 #include "ReconstructionDataFormats/Track.h"
 
 #include "TLorentzVector.h"
 
 #include <vector>
 
@@ -33,8 +33,10 @@
   float mE;
   float mVx, mVy, mVz;
   float mPx, mPy, mPz;
+  bool mIsAlive;
 
   // Setters
+  void setIsAlive(bool isAlive) { mIsAlive = isAlive; }
   void setPDG(int pdg) { mPdgCode = pdg; }
   void setVxVyVz(float vx, float vy, float vz)
   {
@@ -52,6 +54,7 @@
 
   // Getters
   int pdgCode() const { return mPdgCode; }
+  int isAlive() const { return mIsAlive; }
   float vx() const { return mVx; }
   float vy() const { return mVy; }
   float vz() const { return mVz; }
@@ -66,8 +69,8 @@
 /// \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<double> productionVertex,
                                     o2::track::TrackParCov& o2track);
 
@@ -78,8 +81,8 @@
 /// \param o2track the address of the resulting TrackParCov
 /// \param pdg the pdg service
 template <typename PdgService>
 void convertTLorentzVectorToO2Track(int pdgCode,
                                     TLorentzVector particle,
                                     std::vector<double> productionVertex,
                                     o2::track::TrackParCov& o2track,
                                     const PdgService& pdg)
@@ -89,7 +92,7 @@
   if (pdgInfo != nullptr) {
     charge = pdgInfo->Charge() / 3;
   }
   convertTLorentzVectorToO2Track(charge, particle, productionVertex, o2track);
 }
 
 /// Function to convert a OTFParticle into a perfect Track
@@ -97,7 +100,9 @@
 /// \param o2track the address of the resulting TrackParCov
 /// \param pdg the pdg service
 template <typename PdgService>
-void convertOTFParticleToO2Track(const OTFParticle& particle, o2::track::TrackParCov& o2track, const PdgService& pdg)
+void convertOTFParticleToO2Track(const OTFParticle& particle,
+                                 o2::track::TrackParCov& o2track,
+                                 const PdgService& pdg)
 {
   static TLorentzVector tlv;
   tlv.SetPxPyPzE(particle.px(), particle.py(), particle.pz(), particle.e());

ALICE3/DataModel/OTFMCParticle.h

@@ -32,6 +32,9 @@ 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(IsAlive, isAlive, bool);
+DECLARE_SOA_COLUMN(IsPrimary, isPrimary, bool);
+
 } // namespace otfmcparticle
 
 DECLARE_SOA_TABLE_FULL(McParticlesWithDau, "McParticlesWithDau", "AOD", "MCPARTICLEWITHDAU",
@@ -56,6 +59,8 @@ DECLARE_SOA_TABLE_FULL(McParticlesWithDau, "McParticlesWithDau", "AOD", "MCPARTI
                        otfmcparticle::Pt,
                        otfmcparticle::P,
                        otfmcparticle::Y,
+                       otfmcparticle::IsAlive,
+                       otfmcparticle::IsPrimary,
                        mcparticle::PVector<mcparticle::Px, mcparticle::Py, mcparticle::Pz>,
                        mcparticle::ProducedByGenerator<mcparticle::Flags>,
                        mcparticle::FromBackgroundEvent<mcparticle::Flags>,

ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx

@@ -51,7 +51,7 @@ using namespace o2::framework;
 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 constexpr float Tolerance = 1e-7f;
 static const std::vector<std::string> parameterNames{"enable"};
 static const std::vector<std::string> particleNames{"K0s",
                                                     "Lambda",
@@ -102,6 +102,12 @@ struct OnTheFlyDecayer {
     histos.add("AntiXi/hGenAntiXi", "hGenAntiXi", kTH1D, {axisPt});
     histos.add("Omega/hGenOmega", "hGenOmega", kTH1D, {axisPt});
     histos.add("AntiOmega/hGenAntiOmega", "hGenAntiOmega", kTH1D, {axisPt});
+
+    histos.add("GeneratedElectron/hGenEl", "hGenEl", kTH1D, {axisPt});
+    histos.add("GeneratedMuon/hGenMu", "hGenMu", kTH1D, {axisPt});
+    histos.add("GeneratedPion/hGenPi", "hGenPi", kTH1D, {axisPt});
+    histos.add("GeneratedKaon/hGenKa", "hGenKa", kTH1D, {axisPt});
+    histos.add("GeneratedProton/hGenPr", "hGenPr", kTH1D, {axisPt});
   }
 
   bool canDecay(const int pdgCode)
@@ -123,40 +129,60 @@ struct OnTheFlyDecayer {
         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::upgrade::OTFParticle& dau = decayDaughters[idau];
           o2::track::TrackParCov dauTrack;
           o2::upgrade::convertOTFParticleToO2Track(dau, dauTrack, pdgDB);
           if (canDecay(dau.pdgCode())) {
+            dau.setIsAlive(false);
             std::vector<o2::upgrade::OTFParticle> cascadingDaughers = decayer.decayParticle(pdgDB, dauTrack, dau.pdgCode());
             for (const auto& daudau : cascadingDaughers) {
               decayDaughters.push_back(daudau);
               if (MaxNestedDecays < ++nDecays) {
                 LOG(error) << "Seemingly stuck trying to perpetually decay products from pdg: " << particle.pdgCode();
               }
             }
+          } else {
+            dau.setIsAlive(true);
           }
         }
 
         if (decayDaughters.empty()) {
           LOG(error) << "Attempted to decay " << particle.pdgCode() << " but resulting vector of daugthers were empty";
           continue;
         }
-      }
 
-      if (particle.pdgCode() == kK0Short) {
-        histos.fill(HIST("K0S/hGenK0S"), particle.pt());
-      } else if (particle.pdgCode() == kLambda0) {
-        histos.fill(HIST("Lambda/hGenLambda"), particle.pt());
-      } else if (particle.pdgCode() == kLambda0Bar) {
-        histos.fill(HIST("AntiLambda/hGenAntiLambda"), particle.pt());
-      } else if (particle.pdgCode() == kXiMinus) {
-        histos.fill(HIST("Xi/hGenXi"), particle.pt());
-      } else if (particle.pdgCode() == kXiPlusBar) {
-        histos.fill(HIST("AntiXi/hGenAntiXi"), particle.pt());
-      } else if (particle.pdgCode() == kOmegaMinus) {
-        histos.fill(HIST("Omega/hGenOmega"), particle.pt());
-      } else if (particle.pdgCode() == kOmegaPlusBar) {
-        histos.fill(HIST("AntiOmega/hGenAntiOmega"), particle.pt());
+        switch (particle.pdgCode()) {
+          case kK0Short:
+            histos.fill(HIST("K0S/hGenK0S"), particle.pt());
+            break;
+
+          case kLambda0:
+            histos.fill(HIST("Lambda/hGenLambda"), particle.pt());
+            break;
+
+          case kLambda0Bar:
+            histos.fill(HIST("AntiLambda/hGenAntiLambda"), particle.pt());
+            break;
+
+          case kXiMinus:
+            histos.fill(HIST("Xi/hGenXi"), particle.pt());
+            break;
+
+          case kXiPlusBar:
+            histos.fill(HIST("AntiXi/hGenAntiXi"), particle.pt());
+            break;
+
+          case kOmegaMinus:
+            histos.fill(HIST("Omega/hGenOmega"), particle.pt());
+            break;
+
+          case kOmegaPlusBar:
+            histos.fill(HIST("AntiOmega/hGenAntiOmega"), particle.pt());
+            break;
+
+          default:
+            break;
+        }
       }
 
       int daughtersIdSlice[2];
@@ -178,13 +204,13 @@ struct OnTheFlyDecayer {
       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()) < VerySmall) {
+      if ((p - particle.pz()) < Tolerance) {
         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()) < VerySmall) {
+      if ((particle.e() - particle.pz()) < Tolerance) {
         y = (particle.pz() < 0.0f) ? -100.0f : 100.0f;
       } else {
         y = 0.5f * std::log((particle.e() + particle.pz()) / (particle.e() - particle.pz()));
@@ -196,7 +222,7 @@ struct OnTheFlyDecayer {
                               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);
+                              phi, eta, pt, p, y, !canDecay(particle.pdgCode()), true);
     }
 
     int daughtersIdSlice[2] = {-1, -1};
@@ -217,26 +243,51 @@ 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()) < VerySmall) {
+        if ((p - dau.pz()) < Tolerance) {
           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()) < VerySmall) {
+        if ((dau.e() - dau.pz()) < Tolerance) {
           y = (dau.pz() < 0.0f) ? -100.0f : 100.0f;
         } else {
           y = 0.5f * std::log((dau.e() + dau.pz()) / (dau.e() - dau.pz()));
         }
 
+        switch (dau.pdgCode()) {
+          case kElectron:
+            histos.fill(HIST("GeneratedElectron/hGenEl"), pt);
+            break;
+
+          case kMuonMinus:
+            histos.fill(HIST("GeneratedMuon/hGenMu"), pt);
+            break;
+
+          case kPiPlus:
+            histos.fill(HIST("GeneratedPion/hGenPi"), pt);
+            break;
+
+          case kKPlus:
+            histos.fill(HIST("GeneratedKaon/hGenKa"), pt);
+            break;
+
+          case kProton:
+            histos.fill(HIST("GeneratedProton/hGenPr"), pt);
+            break;
+
+          default:
+            break;
+        }
+
         // 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);
+                                phi, eta, pt, p, y, dau.isAlive(), false);
       }
     }
   }