Merge remote-tracking branch 'upstream/master'

Author: sangwoo184

ALICE3/DataModel/OTFMCParticle.h

@@ -37,7 +37,7 @@ DECLARE_SOA_COLUMN(IsPrimary, isPrimary, bool);
 
 } // namespace otfmcparticle
 
-DECLARE_SOA_TABLE_FULL(McParticlesWithDau, "McParticlesWithDau", "AOD", "MCPARTICLEWITHDAU",
+DECLARE_SOA_TABLE_FULL(McPartsWithDau, "McPartsWithDau", "AOD", "MCPARTSWITHDAU",
                        o2::soa::Index<>,
                        mcparticle::McCollisionId,
                        mcparticle::PdgCode,
@@ -69,7 +69,7 @@ DECLARE_SOA_TABLE_FULL(McParticlesWithDau, "McParticlesWithDau", "AOD", "MCPARTI
                        mcparticle::GetProcess<mcparticle::Flags, mcparticle::StatusCode>,
                        mcparticle::IsPhysicalPrimary<mcparticle::Flags>);
 
-using McParticleWithDau = McParticlesWithDau::iterator;
+using McPartWithDau = McPartsWithDau::iterator;
 
 } // namespace o2::aod
 

ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx

@@ -19,23 +19,13 @@
 #include "ALICE3/Core/TrackUtilities.h"
 #include "ALICE3/DataModel/OTFMCParticle.h"
 
-#include <CCDB/BasicCCDBManager.h>
 #include <CommonConstants/MathConstants.h>
 #include <CommonConstants/PhysicsConstants.h>
-#include <DCAFitter/DCAFitterN.h>
-#include <DataFormatsParameters/GRPMagField.h>
-#include <DetectorsBase/Propagator.h>
-#include <DetectorsVertexing/PVertexer.h>
-#include <DetectorsVertexing/PVertexerHelpers.h>
-#include <Field/MagneticField.h>
 #include <Framework/AnalysisDataModel.h>
 #include <Framework/AnalysisTask.h>
 #include <Framework/HistogramRegistry.h>
 #include <Framework/O2DatabasePDGPlugin.h>
-#include <Framework/StaticFor.h>
 #include <Framework/runDataProcessing.h>
-#include <ReconstructionDataFormats/DCA.h>
-#include <SimulationDataFormat/InteractionSampler.h>
 
 #include <TPDGCode.h>
 
@@ -69,14 +59,15 @@ static const std::vector<int> pdgCodes{kK0Short,
                                        kOmegaPlusBar};
 
 struct OnTheFlyDecayer {
-  Produces<aod::McParticlesWithDau> tableMcParticlesWithDau;
+  Produces<aod::McPartsWithDau> tableMcParticlesWithDau;
 
   o2::upgrade::Decayer decayer;
   Service<o2::framework::O2DatabasePDG> pdgDB;
   std::map<int, std::vector<o2::upgrade::OTFParticle>> mDecayDaughters;
 
   Configurable<int> seed{"seed", 0, "Set seed for particle decayer"};
   Configurable<float> magneticField{"magneticField", 20., "Magnetic field (kG)"};
+  Configurable<float> maxEta{"maxEta", 2.5, "Only decay particles that appear within selected eta range"};
   Configurable<LabeledArray<int>> enabledDecays{"enabledDecays",
                                                 {DefaultParameters[0], NumDecays, NumParameters, particleNames, parameterNames},
                                                 "Enable option for particle to be decayed: 0 - no, 1 - yes"};
@@ -122,6 +113,15 @@ struct OnTheFlyDecayer {
                                                      y(y),
                                                      isAlive(isAlive),
                                                      isPrimary(isPrimary) {}
+
+    bool hasNaN() const
+    {
+      return std::isnan(px) || std::isnan(py) || std::isnan(pz) || std::isnan(e) ||
+             std::isnan(vx) || std::isnan(vy) || std::isnan(vz) || std::isnan(vt) ||
+             std::isnan(phi) || std::isnan(eta) || std::isnan(pt) || std::isnan(p) ||
+             std::isnan(y) || std::isnan(weight);
+    }
+
     int collisionId;
     int pdgCode;
     int statusCode;
@@ -148,6 +148,10 @@ struct OnTheFlyDecayer {
       }
     }
 
+    auto hNaNBookkeeping = histos.add<TH1>("hNaNBookkeeping", "hNaNBookkeeping", kTH1D, {{2, -0.5, 1.5}});
+    hNaNBookkeeping->GetXaxis()->SetBinLabel(1, "OK");
+    hNaNBookkeeping->GetXaxis()->SetBinLabel(2, "NaN");
+
     histos.add("K0S/hGenK0S", "hGenK0S", kTH1D, {axisPt});
     histos.add("Lambda/hGenLambda", "hGenLambda", kTH1D, {axisPt});
     histos.add("AntiLambda/hGenAntiLambda", "hGenAntiLambda", kTH1D, {axisPt});
@@ -175,30 +179,29 @@ struct OnTheFlyDecayer {
     mcParticlesAlice3.clear();
     u_int64_t nStoredDaughters = 0;
     for (int index{0}; index < static_cast<int>(mcParticles.size()); ++index) {
-      const auto& particle = mcParticles.iteratorAt(index);
-      std::vector<o2::upgrade::OTFParticle> decayDaughters;
-      static constexpr int MaxNestedDecays = 10;
-      int nDecays = 0;
-      if (canDecay(particle.pdgCode())) {
+      const auto& particle = mcParticles.rawIteratorAt(index);
+      std::vector<o2::upgrade::OTFParticle> decayDaughters, decayStack;
+      if (canDecay(particle.pdgCode()) && std::abs(particle.eta()) < maxEta) {
         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) {
-          o2::upgrade::OTFParticle dau = decayDaughters[idau];
+        decayStack = decayer.decayParticle(pdgDB, o2track, particle.pdgCode());
+        while (!decayStack.empty()) {
+          o2::upgrade::OTFParticle otfParticle = decayStack.back();
+          decayStack.pop_back();
+
+          const bool stable = !canDecay(otfParticle.pdgCode());
+          otfParticle.setIsAlive(stable);
+          decayDaughters.push_back(otfParticle);
+
+          if (stable) {
+            continue;
+          }
+
           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 (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();
-              }
-            }
-          } else {
-            dau.setIsAlive(true);
+          o2::upgrade::convertOTFParticleToO2Track(otfParticle, dauTrack, pdgDB);
+          std::vector<o2::upgrade::OTFParticle> daughters = decayer.decayParticle(pdgDB, dauTrack, otfParticle.pdgCode());
+          for (o2::upgrade::OTFParticle dau : daughters) {
+            decayStack.push_back(dau);
           }
         }
 
@@ -338,7 +341,7 @@ struct OnTheFlyDecayer {
         // TODO: Particle status code
         // TODO: Expression columns
         // TODO: vt
-        auto mother = mcParticles.iteratorAt(index);
+        auto mother = mcParticles.rawIteratorAt(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(),
@@ -348,8 +351,13 @@ struct OnTheFlyDecayer {
     }
 
     for (const auto& particle : mcParticlesAlice3) {
-      std::span<const int> motherSpan(particle.mothersIds, 2);
+      if (particle.hasNaN()) {
+        histos.fill(HIST("hNaNBookkeeping"), 1);
+        continue;
+      }
 
+      histos.fill(HIST("hNaNBookkeeping"), 0);
+      std::span<const int> 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,

ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

@@ -51,6 +51,7 @@
 #include <Framework/StaticFor.h>
 #include <Framework/runDataProcessing.h>
 #include <ReconstructionDataFormats/DCA.h>
+#include <ReconstructionDataFormats/PID.h>
 #include <SimulationDataFormat/InteractionSampler.h>
 
 #include <TGenPhaseSpace.h>
@@ -597,6 +598,31 @@ struct OnTheFlyTracker {
     gRandom->SetSeed(seed);
   }
 
+  /// Function to get the internal PID for a given pdgCode
+  /// \param pdgCode pdg code for a common particle (particle handled by the tracker)
+  int pdgCodeToPID(int pdgCode) const
+  {
+    if (std::abs(pdgCode) == PDG_t::kElectron) {
+      return o2::track::PID::Electron;
+    } else if (std::abs(pdgCode) == PDG_t::kMuonMinus) {
+      return o2::track::PID::Muon;
+    } else if (std::abs(pdgCode) == PDG_t::kPiPlus) {
+      return o2::track::PID::Pion;
+    } else if (std::abs(pdgCode) == PDG_t::kKPlus) {
+      return o2::track::PID::Kaon;
+    } else if (std::abs(pdgCode) == PDG_t::kProton) {
+      return o2::track::PID::Proton;
+    } else if (std::abs(pdgCode) == PDG_t::kLambda0) {
+      return o2::track::PID::Lambda;
+    } else if (std::abs(pdgCode) == PDG_t::kXiMinus) {
+      return o2::track::PID::XiMinus;
+    } else if (std::abs(pdgCode) == PDG_t::kOmegaMinus) {
+      return o2::track::PID::OmegaMinus;
+    } else {
+      return o2::track::PID::Pion; // Default trackParCov assumption
+    }
+  }
+
   /// Function to decay the xi
   /// \param particle the particle to decay
   /// \param track track of particle to decay
@@ -705,7 +731,7 @@ struct OnTheFlyTracker {
 
   float dNdEta = 0.f;                                                      // Charged particle multiplicity to use in the efficiency evaluation
   std::pair<float, float> vertexReconstructionEfficiencyCounters = {0, 0}; // {nVerticesWithMoreThan2Contributors, nVerticesReconstructed}
-  void processWithLUTs(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, int const& icfg)
+  void processWithLUTs(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const int icfg)
   {
     vertexReconstructionEfficiencyCounters.first += 1;
     const int lastTrackIndex = tableStoredTracksCov.lastIndex() + 1; // bookkeep the last added track
@@ -878,9 +904,12 @@ struct OnTheFlyTracker {
           getHist(TH1, histPath + "hXiBuilding")->Fill(0.0f);
         }
 
-        o2::upgrade::convertTLorentzVectorToO2Track(kPiMinus, cascadeDecayProducts[0], xiDecayVertex, xiDaughterTrackParCovsPerfect[0], pdgDB);
-        o2::upgrade::convertTLorentzVectorToO2Track(kPiMinus, cascadeDecayProducts[1], laDecayVertex, xiDaughterTrackParCovsPerfect[1], pdgDB);
-        o2::upgrade::convertTLorentzVectorToO2Track(kProton, cascadeDecayProducts[2], laDecayVertex, xiDaughterTrackParCovsPerfect[2], pdgDB);
+        o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[0], xiDecayVertex, xiDaughterTrackParCovsPerfect[0], pdgDB);
+        xiDaughterTrackParCovsPerfect[0].setPID(pdgCodeToPID(PDG_t::kPiMinus));
+        o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[1], laDecayVertex, xiDaughterTrackParCovsPerfect[1], pdgDB);
+        xiDaughterTrackParCovsPerfect[1].setPID(pdgCodeToPID(PDG_t::kPiMinus));
+        o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kProton, cascadeDecayProducts[2], laDecayVertex, xiDaughterTrackParCovsPerfect[2], pdgDB);
+        xiDaughterTrackParCovsPerfect[2].setPID(pdgCodeToPID(PDG_t::kProton));
 
         for (int i = 0; i < kCascProngs; i++) {
           isReco[i] = false;
@@ -1001,7 +1030,7 @@ struct OnTheFlyTracker {
               {posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]},
               covV, 0, true);
             v0Track.setAbsCharge(0);
-            v0Track.setPID(o2::track::PID::Lambda);
+            v0Track.setPID(pdgCodeToPID(PDG_t::kLambda0));
 
             // dca fitter step
             nCand = 0;
@@ -1040,8 +1069,8 @@ struct OnTheFlyTracker {
 
               // initialize cascade track
               o2::track::TrackParCov cascadeTrack = fitter.createParentTrackParCov();
-              cascadeTrack.setAbsCharge(-1);                // may require more adjustments
-              cascadeTrack.setPID(o2::track::PID::XiMinus); // FIXME: not OK for omegas
+              cascadeTrack.setAbsCharge(-1);                      // may require more adjustments
+              cascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas
 
               thisCascade.cascradiusMC = xiDecayRadius2D;
               thisCascade.pt = cascadeTrack.getPt();
@@ -1131,6 +1160,7 @@ struct OnTheFlyTracker {
           o2::upgrade::convertMCParticleToO2Track(mcParticle, prefectCascadeTrack, pdgDB);
 
           // back track is already smeared
+          prefectCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas
           int nCascHits = fastTracker[icfg]->FastTrack(prefectCascadeTrack, trackedCasc, dNdEta);
           reconstructedCascade = (fastTrackerSettings.minSiliconHitsForKinkReco < nCascHits) ? false : true;
 
@@ -1201,6 +1231,7 @@ struct OnTheFlyTracker {
                                                       std::array{pV0[0], pV0[1], pV0[2]}},
                                            std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda});
 
+            newCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas
             tracksAlice3.push_back(TrackAlice3{newCascadeTrack, mcParticle.globalIndex(), time, timeResolutionUs, false, false, 1, thisCascade.foundClusters});
 
             // add this cascade to vector (will fill cursor later with collision ID)
@@ -1436,6 +1467,7 @@ struct OnTheFlyTracker {
       } else if (fastPrimaryTrackerSettings.fastTrackPrimaries) {
         o2::track::TrackParCov o2Track;
         o2::upgrade::convertMCParticleToO2Track(mcParticle, o2Track, pdgDB);
+        o2Track.setPID(pdgCodeToPID(mcParticle.pdgCode()));
         const int nHits = fastTracker[icfg]->FastTrack(o2Track, trackParCov, dNdEta);
         if (nHits < fastPrimaryTrackerSettings.minSiliconHits) {
           reconstructed = false;
@@ -1736,7 +1768,7 @@ struct OnTheFlyTracker {
     }
   }
 
-  void processConfigurationDev(aod::McCollision const& mcCollision, aod::McParticlesWithDau const& mcParticles, int const& icfg)
+  void processConfigurationDev(aod::McCollision const& mcCollision, aod::McPartsWithDau const& mcParticles, const int icfg)
   {
     // const int lastTrackIndex = tableStoredTracksCov.lastIndex() + 1; // bookkeep the last added track
     const std::string histPath = "Configuration_" + std::to_string(icfg) + "/";
@@ -1838,6 +1870,7 @@ struct OnTheFlyTracker {
       } else if (enableSecondarySmearing) {
         o2::track::TrackParCov perfectTrackParCov;
         o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectTrackParCov, pdgDB);
+        perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
         const int nHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta);
         if (nHits < fastTrackerSettings.minSiliconHits) {
           reconstructed = false;
@@ -1991,7 +2024,7 @@ struct OnTheFlyTracker {
     }
   }
 
-  void processDecayer(aod::McCollision const& mcCollision, aod::McParticlesWithDau const& mcParticles)
+  void processDecayer(aod::McCollision const& mcCollision, aod::McPartsWithDau const& mcParticles)
   {
     for (size_t icfg = 0; icfg < mSmearer.size(); ++icfg) {
       processConfigurationDev(mcCollision, mcParticles, static_cast<int>(icfg));

ALICE3/Tasks/CMakeLists.txt

@@ -95,6 +95,6 @@ o2physics_add_dpl_workflow(alice3-dq-efficiency
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(alice3-decayer-qa
-                    SOURCES alice3DecayerQA.cxx
+                    SOURCES alice3DecayerQa.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)