[ALICE3] Propagate more information from otf decayer to otf tracker (#17659)

Author: jesgum

ALICE3/Core/OTFParticle.h

@@ -78,6 +78,7 @@ class OTFParticle
   void setIndicesDaughter(const int start, const int stop) { mIndicesDaughter = {start, stop}; }
   void setProductionTime(const float vt) { mVt = vt; }
   void setFlags(uint8_t flag) { mFlag = flag; }
+  void setDecayRadius(const float decayRadius) { mDecayRadius = decayRadius; }
   void setVxVyVz(const float vx, const float vy, const float vz)
   {
     mVx = vx;
@@ -123,6 +124,7 @@ class OTFParticle
   float pz() const { return mPz; }
   float e() const { return mE; }
   float radius() const { return std::hypot(mVx, mVy); }
+  float decayRadius() const { return mDecayRadius; }
   float r() const { return radius(); }
   float pt() const { return std::hypot(mPx, mPy); }
   float p() const { return std::hypot(mPx, mPy, mPz); }
@@ -185,6 +187,7 @@ class OTFParticle
   int mCollisionId{-1};
   float mVx{}, mVy{}, mVz{}, mVt{};
   float mPx{}, mPy{}, mPz{}, mE{};
+  float mDecayRadius{-1};
 
   int mStatusCode{};
   uint8_t mFlag{};

ALICE3/DataModel/tracksAlice3.h

@@ -50,16 +50,20 @@ using TrackExtraA3 = TracksExtraA3::iterator;
 
 namespace mcparticle_alice3
 {
-DECLARE_SOA_COLUMN(NHits, nHits, int);     //! number of silicon hits
-DECLARE_SOA_COLUMN(Charge, charge, float); //! particle charge
+DECLARE_SOA_COLUMN(NHits, nHits, int);                  //! number of silicon hits
+DECLARE_SOA_COLUMN(Charge, charge, float);              //! particle charge
+DECLARE_SOA_COLUMN(DecayRadius, decayRadius, float);    //! Radius for decayed particle produced by decayer (-1 if not decayed)
 DECLARE_SOA_BITMAP_COLUMN(DecayerBits, decayerBits, 8); //! Bit mask for particle produced by the OTF decayer
 } // namespace mcparticle_alice3
 DECLARE_SOA_TABLE(MCParticlesExtraA3, "AOD", "MCParticlesExtraA3",
                   mcparticle_alice3::NHits,
                   mcparticle_alice3::Charge);
 using MCParticleExtraA3 = MCParticlesExtraA3::iterator;
 
-DECLARE_SOA_TABLE(OTFDecayerBits, "AOD", "OTFDecayerBits", mcparticle_alice3::DecayerBits);
+DECLARE_SOA_TABLE(OTFParticleExtras, "AOD", "OTFPARTICLEEXTRA",
+                  mcparticle_alice3::DecayerBits,
+                  mcparticle_alice3::DecayRadius);
+using OTFParticleExtra = OTFParticleExtras::iterator;
 
 } // namespace o2::aod
 

ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx

@@ -88,7 +88,7 @@ O2ORIGIN("TMP");
 struct OnTheFlyDecayer {
   Produces<aod::McCollisions_001> tableMcCollisions;
   Produces<aod::StoredMcParticles_001> tableMcParticles;
-  Produces<aod::OTFDecayerBits> tableOTFDecayerBits;
+  Produces<aod::OTFParticleExtras> tableOTFParticleExtras;
 
   o2::upgrade::Decayer decayer;
   Service<o2::framework::O2DatabasePDG> pdgDB{};
@@ -172,6 +172,7 @@ struct OnTheFlyDecayer {
         trackLength = o2::upgrade::computeTrackLength(o2track, decayRadius, magneticField);
       }
 
+      particle.setDecayRadius(std::hypot(decayer.getSecondaryVertexX(), decayer.getSecondaryVertexY()));
       const float trackTimeNS = trackLength / trackVelocity * PicoToNano;
       particle.setIndicesDaughter(particlesInDataframe - indexOffset + allParticles.size(), particlesInDataframe - indexOffset + allParticles.size() + (decayStack.size() - 1));
       for (auto& daughter : decayStack) {
@@ -220,7 +221,7 @@ struct OnTheFlyDecayer {
         histos.fill(HIST("hNaNBookkeeping"), 0);
       }
 
-      tableOTFDecayerBits(otfParticle.getBitsValue());
+      tableOTFParticleExtras(otfParticle.getBitsValue(), otfParticle.decayRadius());
       tableMcParticles(tableMcCollisions.lastIndex(), otfParticle.pdgCode(), otfParticle.statusCode(), otfParticle.flags(),
                        otfParticle.getMotherSpan(), otfParticle.getDaughters().data(), otfParticle.weight(),
                        otfParticle.px(), otfParticle.py(), otfParticle.pz(), otfParticle.e(),

ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

@@ -619,6 +619,14 @@ struct OnTheFlyTracker {
         insertHist(histPath + "h2dSecondaryKaPtRes", "h2dSecondaryKaPtRes;Gen p_{T};#Delta p_{T} / Reco p_{T}", kTH2D, {{axes.axisMomentum, axes.axisPtRes}});
         insertHist(histPath + "h2dPrimaryPrPtRes", "h2dPrimaryPrPtRes;Gen p_{T};#Delta p_{T} / Reco p_{T}", kTH2D, {{axes.axisMomentum, axes.axisPtRes}});
         insertHist(histPath + "h2dSecondaryPrPtRes", "h2dSecondaryPrPtRes;Gen p_{T};#Delta p_{T} / Reco p_{T}", kTH2D, {{axes.axisMomentum, axes.axisPtRes}});
+
+        if (fastPrimaryTrackerSettings.fastTrackShortLivedParticles) {
+          insertHist(histPath + "h2dGenShortLivedParticleRadius", "h2dGenShortLivedParticleRadius;Radius (cm);Momentum p_{T}", kTH2D, {{axes.axisDecayRadius, axes.axisMomentum}});
+          insertHist(histPath + "h2dRecShortLivedParticleRadius", "h2dRecShortLivedParticleRadius;Radius (cm);Momentum p_{T}", kTH2D, {{axes.axisDecayRadius, axes.axisMomentum}});
+          insertHist(histPath + "h2dGenRadiusIniVsDecay", "h2dGenRadiusIniVsDecay;Radius (cm);Radius (cm)", kTH2D, {{axes.axisDecayRadius, axes.axisDecayRadius}});
+          insertHist(histPath + "h2dRecRadiusIniVsDecay", "h2dRecRadiusIniVsDecay;Radius (cm);Radius (cm)", kTH2D, {{axes.axisDecayRadius, axes.axisDecayRadius}});
+          insertHist(histPath + "h2dDecayRadiusVsNhits", "h2dDecayRadiusVsNhits;Radius (cm);Nhits", kTH2D, {{axes.axisDecayRadius, {20, 0.5, 20}}});
+        }
       }
 
     } // end config loop
@@ -2151,11 +2159,17 @@ struct OnTheFlyTracker {
         o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectTrackParCov, pdgDB);
         perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
         computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov);
-        nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta);
+        nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta, mcParticle.decayRadius());
+        getHist<TH2>(histPath + "h2dGenShortLivedParticleRadius")->Fill(mcParticle.decayRadius(), perfectTrackParCov.getPt());
+        getHist<TH2>(histPath + "h2dGenRadiusIniVsDecay")->Fill(std::hypot(perfectTrackParCov.getX(), perfectTrackParCov.getY()), mcParticle.decayRadius());
         if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) {
           reconstructed = false;
         } else {
           reconstructed = true;
+          getHist<TH2>(histPath + "h2dRecShortLivedParticleRadius")->Fill(mcParticle.decayRadius(), perfectTrackParCov.getPt());
+          getHist<TH2>(histPath + "h2dRecRadiusIniVsDecay")->Fill(std::hypot(perfectTrackParCov.getX(), perfectTrackParCov.getY()), mcParticle.decayRadius());
+          getHist<TH2>(histPath + "h2dDecayRadiusVsNhits")->Fill(mcParticle.decayRadius(), nTrkHits);
+          LOG(info) << mcParticle.decayRadius();
         }
       } else if (enableSecondarySmearing && isSecondary) {
         o2::track::TrackParCov perfectTrackParCov;
@@ -2259,7 +2273,7 @@ struct OnTheFlyTracker {
     fillTracksInfo(ghostTracksAlice3, primaryVertex, icfg);
   }
 
-  void processDecayer(aod::McCollision const& mcCollision, soa::Join<aod::McParticles, aod::OTFDecayerBits> const& mcParticles)
+  void processDecayer(aod::McCollision const& mcCollision, soa::Join<aod::McParticles, aod::OTFParticleExtras> const& mcParticles)
   {
     for (size_t icfg = 0; icfg < mSmearer.size(); ++icfg) {
       processConfigurationDev(mcCollision, mcParticles, static_cast<int>(icfg));