diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h
index 2280261dd21..51296ffe67d 100644
--- a/ALICE3/Core/Decayer.h
+++ b/ALICE3/Core/Decayer.h
@@ -19,6 +19,7 @@
 #ifndef ALICE3_CORE_DECAYER_H_
 #define ALICE3_CORE_DECAYER_H_
 
+#include "ALICE3/Core/OTFParticle.h"
 #include "ALICE3/Core/TrackUtilities.h"
 
 #include <CommonConstants/PhysicsConstants.h>
@@ -60,13 +61,12 @@ class Decayer
     const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm
     const double betaGamma = particle.p() / mass;
     const double rxyz = -betaGamma * ctau * std::log(1 - u);
-    double vx, vy, vz;
     double px, py, e;
 
     if (!charge) {
-      vx = particle.vx() + rxyz * (particle.px() / particle.p());
-      vy = particle.vy() + rxyz * (particle.py() / particle.p());
-      vz = particle.vz() + rxyz * (particle.pz() / particle.p());
+      mVx = particle.vx() + rxyz * (particle.px() / particle.p());
+      mVy = particle.vy() + rxyz * (particle.py() / particle.p());
+      mVz = particle.vz() + rxyz * (particle.pz() / particle.p());
       px = particle.px();
       py = particle.py();
     } else {
@@ -74,14 +74,14 @@ class Decayer
       o2::math_utils::CircleXYf_t circle;
       o2::upgrade::convertOTFParticleToO2Track(particle, track, pdgDB);
 
-      float sna, csa;
+      float sna{}, csa{};
       track.getCircleParams(mBz, circle, sna, csa);
       const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl());
       const double theta = rxy / circle.rC;
 
-      vx = ((particle.vx() - circle.xC) * std::cos(theta) - (particle.vy() - circle.yC) * std::sin(theta)) + circle.xC;
-      vy = ((particle.vy() - circle.yC) * std::cos(theta) + (particle.vx() - circle.xC) * std::sin(theta)) + circle.yC;
-      vz = particle.vz() + rxyz * (particle.pz() / track.getP());
+      mVx = ((particle.vx() - circle.xC) * std::cos(theta) - (particle.vy() - circle.yC) * std::sin(theta)) + circle.xC;
+      mVy = ((particle.vy() - circle.yC) * std::cos(theta) + (particle.vx() - circle.xC) * std::sin(theta)) + circle.yC;
+      mVz = particle.vz() + rxyz * (particle.pz() / track.getP());
 
       px = particle.px() * std::cos(theta) - particle.py() * std::sin(theta);
       py = particle.py() * std::cos(theta) + particle.px() * std::sin(theta);
@@ -124,7 +124,7 @@ class Decayer
       o2::upgrade::OTFParticle particle;
       TLorentzVector dau = *decay.GetDecay(i);
       particle.setPDG(pdgCodesDaughters[i]);
-      particle.setVxVyVz(vx, vy, vz);
+      particle.setVxVyVz(mVx, mVy, mVz);
       particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E());
       decayProducts.push_back(particle);
     }
@@ -132,17 +132,24 @@ class Decayer
     return decayProducts;
   }
 
+  // Setters
+  void setBField(const double b) { mBz = b; }
   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
+  float getSecondaryVertexX() const { return static_cast<float>(mVx); }
+  float getSecondaryVertexY() const { return static_cast<float>(mVy); }
+  float getSecondaryVertexZ() const { return static_cast<float>(mVz); }
+  float getDecayRadius() const { return static_cast<float>(std::hypot(mVx, mVy)); }
 
  private:
-  TRandom3 mRand3;
-  double mBz;
+  double mBz{20.}; // kG
+  double mVx{-1.}, mVy{-1.}, mVz{-1.};
+  TRandom3 mRand3{};
 };
 
 } // namespace upgrade
diff --git a/ALICE3/Core/OTFParticle.h b/ALICE3/Core/OTFParticle.h
new file mode 100644
index 00000000000..b6ca1f246c3
--- /dev/null
+++ b/ALICE3/Core/OTFParticle.h
@@ -0,0 +1,165 @@
+// 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 OTFParticle.h
+/// \brief Basic class to hold information regarding a mc particle to be used in fast simulation
+/// \author Jesper Karlsson Gumprecht <jesper.gumprecht@cern.ch>
+///
+
+#ifndef ALICE3_CORE_OTFPARTICLE_H_
+#define ALICE3_CORE_OTFPARTICLE_H_
+
+#include <CommonConstants/MathConstants.h>
+
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <span>
+
+namespace o2::upgrade
+{
+
+class OTFParticle
+{
+ public:
+  OTFParticle() = default;
+
+  template <typename TParticle>
+  explicit OTFParticle(const TParticle& particle)
+  {
+    mPdgCode = particle.pdgCode();
+    mGlobalIndex = particle.globalIndex();
+    mCollisionId = particle.mcCollisionId();
+    mPx = particle.px();
+    mPy = particle.py();
+    mPz = particle.pz();
+    mE = particle.e();
+    mVx = particle.vx();
+    mVy = particle.vy();
+    mVz = particle.vz();
+    mVt = particle.vt();
+    mIsFromMcParticles = true;
+    if (particle.has_mothers()) {
+      mIndicesMother = {particle.mothersIds().front(), particle.mothersIds().back()};
+    }
+  }
+
+  // Setters
+  void setIsAlive(const bool isAlive) { mIsAlive = isAlive; }
+  void setIsPrimary(const bool isPrimary) { mIsPrimary = isPrimary; }
+  void setCollisionId(const int collisionId) { mCollisionId = collisionId; }
+  void setPDG(const int pdg) { mPdgCode = pdg; }
+  void setIndicesMother(const int start, const int stop) { mIndicesMother = {start, stop}; }
+  void setIndicesDaughter(const int start, const int stop) { mIndicesDaughter = {start, stop}; }
+  void setProductionTime(const float vt) { mVt = vt; }
+  void setVxVyVz(const float vx, const float vy, const float vz)
+  {
+    mVx = vx;
+    mVy = vy;
+    mVz = vz;
+  }
+  void setPxPyPzE(const float px, const float py, const float pz, const float e)
+  {
+    mPx = px;
+    mPy = py;
+    mPz = pz;
+    mE = e;
+  }
+
+  // Getters
+  int pdgCode() const { return mPdgCode; }
+  int globalIndex() const { return mGlobalIndex; }
+  int collisionId() const { return mCollisionId; }
+  bool isAlive() const { return mIsAlive; }
+  bool isPrimary() const { return mIsPrimary; }
+  bool isFromMcParticles() const { return mIsFromMcParticles; }
+  float weight() const
+  {
+    static constexpr float Weight = 1.f;
+    return Weight;
+  }
+  uint8_t flags() const
+  {
+    static constexpr uint8_t Flags = 1;
+    return Flags; // todo
+  }
+  int statusCode() const
+  {
+    static constexpr int StatusCode = 1;
+    return StatusCode; // todo
+  }
+  float vx() const { return mVx; }
+  float vy() const { return mVy; }
+  float vz() const { return mVz; }
+  float vt() const { return mVt; }
+  float px() const { return mPx; }
+  float py() const { return mPy; }
+  float pz() const { return mPz; }
+  float e() const { return mE; }
+  float radius() const { return std::hypot(mVx, mVy); }
+  float r() const { return radius(); }
+  float pt() const { return std::hypot(mPx, mPy); }
+  float p() const { return std::hypot(mPx, mPy, mPz); }
+  float phi() const { return o2::constants::math::PI + std::atan2(-1.0f * py(), -1.0f * px()); }
+  float eta() const
+  {
+    // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943
+    static constexpr float Tolerance = 1e-7f;
+    if ((p() - mPz) < Tolerance) {
+      return (mPz < 0.0f) ? -100.0f : 100.0f;
+    } else {
+      return 0.5f * std::log((p() + mPz) / (p() - mPz));
+    }
+  }
+  float y() const
+  {
+    // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1922
+    static constexpr float Tolerance = 1e-7f;
+    if ((e() - mPz) < Tolerance) {
+      return (mPz < 0.0f) ? -100.0f : 100.0f;
+    } else {
+      return 0.5f * std::log((mE + mPz) / (mE - mPz));
+    }
+  }
+  int getMotherIndexStart() const { return mIndicesMother[0]; }
+  int getMotherIndexStop() const { return mIndicesMother[1]; }
+  int getDaughterIndexStart() const { return mIndicesDaughter[0]; }
+  int getDaughterIndexStop() const { return mIndicesDaughter[1]; }
+  std::array<int, 2> getMothers() const { return mIndicesMother; }
+  std::array<int, 2> getDaughters() const { return mIndicesDaughter; }
+  std::span<const int> getMotherSpan() const { return hasMothers() ? std::span<const int>(mIndicesMother.data(), 2) : std::span<const int>(); }
+
+  // Checks
+  bool hasDaughters() const { return (mIndicesDaughter[0] > 0); }
+  bool hasMothers() const { return (mIndicesMother[0] > 0); }
+  bool hasNaN() const
+  {
+    return std::isnan(mPx) || std::isnan(mPy) || std::isnan(mPz) || std::isnan(mE) ||
+           std::isnan(mVx) || std::isnan(mVy) || std::isnan(mVz);
+  }
+  bool hasIndex() const
+  {
+    return (mGlobalIndex != -1);
+  }
+
+ private:
+  int mPdgCode{}, mGlobalIndex{-1};
+  int mCollisionId{};
+  float mVx{}, mVy{}, mVz{}, mVt{};
+  float mPx{}, mPy{}, mPz{}, mE{};
+  bool mIsAlive{}, mIsFromMcParticles{false};
+  bool mIsPrimary{};
+  std::array<int, 2> mIndicesMother{-1, -1}, mIndicesDaughter{-1, -1};
+};
+
+} // namespace o2::upgrade
+#endif // ALICE3_CORE_OTFPARTICLE_H_
diff --git a/ALICE3/Core/TrackUtilities.cxx b/ALICE3/Core/TrackUtilities.cxx
index 739e1df28ae..4239864e933 100644
--- a/ALICE3/Core/TrackUtilities.cxx
+++ b/ALICE3/Core/TrackUtilities.cxx
@@ -45,3 +45,71 @@ void o2::upgrade::convertTLorentzVectorToO2Track(const int charge,
   // Initialize TrackParCov in-place
   new (&o2track)(o2::track::TrackParCov)(x, particle.Phi(), params, covm);
 }
+
+float o2::upgrade::computeParticleVelocity(float momentum, float mass)
+{
+  const float a = momentum / mass;
+  // uses light speed in cm/ps so output is in those units
+  return o2::constants::physics::LightSpeedCm2PS * a / std::sqrt((1.f + a * a));
+}
+
+float o2::upgrade::computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField)
+{
+  // don't make use of the track parametrization
+  float length = -100;
+
+  o2::math_utils::CircleXYf_t trcCircle;
+  float sna, csa;
+  track.getCircleParams(magneticField, trcCircle, sna, csa);
+
+  // distance between circle centers (one circle is at origin -> easy)
+  const float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC);
+
+  // condition of circles touching - if not satisfied returned length will be -100
+  if (centerDistance < trcCircle.rC + radius && centerDistance > std::fabs(trcCircle.rC - radius)) {
+    length = 0.0f;
+
+    // base radical direction
+    const float ux = trcCircle.xC / centerDistance;
+    const float uy = trcCircle.yC / centerDistance;
+    // calculate perpendicular vector (normalized) for +/- displacement
+    const float vx = -uy;
+    const float vy = +ux;
+    // calculate coordinate for radical line
+    const float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + radius * radius) / (2.0f * centerDistance);
+    // calculate absolute displacement from center-to-center axis
+    const float displace = (0.5f / centerDistance) * std::sqrt(
+                                                       (-centerDistance + trcCircle.rC - radius) *
+                                                       (-centerDistance - trcCircle.rC + radius) *
+                                                       (-centerDistance + trcCircle.rC + radius) *
+                                                       (centerDistance + trcCircle.rC + radius));
+
+    // possible intercept points of track and TOF layer in 2D plane
+    const float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy};
+    const float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy};
+
+    // decide on correct intercept point
+    std::array<float, 3> mom;
+    track.getPxPyPzGlo(mom);
+    const float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1];
+    const float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1];
+
+    // get start point
+    std::array<float, 3> startPoint;
+    track.getXYZGlo(startPoint);
+
+    float cosAngle = -1000, modulus = -1000;
+
+    if (scalarProduct1 > scalarProduct2) {
+      modulus = std::hypot(point1[0] - trcCircle.xC, point1[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
+      cosAngle = (point1[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point1[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC);
+    } else {
+      modulus = std::hypot(point2[0] - trcCircle.xC, point2[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
+      cosAngle = (point2[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point2[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC);
+    }
+    cosAngle /= modulus;
+    length = trcCircle.rC * std::acos(cosAngle);
+    length *= std::sqrt(1.0f + track.getTgl() * track.getTgl());
+  }
+  return length;
+}
diff --git a/ALICE3/Core/TrackUtilities.h b/ALICE3/Core/TrackUtilities.h
index 6e4905672bd..b80ea1b2341 100644
--- a/ALICE3/Core/TrackUtilities.h
+++ b/ALICE3/Core/TrackUtilities.h
@@ -18,6 +18,8 @@
 #ifndef ALICE3_CORE_TRACKUTILITIES_H_
 #define ALICE3_CORE_TRACKUTILITIES_H_
 
+#include "ALICE3/Core/OTFParticle.h"
+
 #include <ReconstructionDataFormats/Track.h>
 
 #include <TLorentzVector.h>
@@ -30,145 +32,6 @@
 namespace o2::upgrade
 {
 
-class OTFParticle
-{
- public:
-  OTFParticle() = default;
-
-  template <typename TParticle>
-  explicit OTFParticle(const TParticle& particle)
-  {
-    mPdgCode = particle.pdgCode();
-    mGlobalIndex = particle.globalIndex();
-    mCollisionId = particle.mcCollisionId();
-    mPx = particle.px();
-    mPy = particle.py();
-    mPz = particle.pz();
-    mE = particle.e();
-    mVx = particle.vx();
-    mVy = particle.vy();
-    mVz = particle.vz();
-    mIsFromMcParticles = true;
-    if (particle.has_mothers()) {
-      mIndicesMother = {particle.mothersIds().front(), particle.mothersIds().back()};
-    }
-  }
-
-  // Setters
-  void setIsAlive(const bool isAlive) { mIsAlive = isAlive; }
-  void setIsPrimary(const bool isPrimary) { mIsPrimary = isPrimary; }
-  void setCollisionId(const int collisionId) { mCollisionId = collisionId; }
-  void setPDG(const int pdg) { mPdgCode = pdg; }
-  void setIndicesMother(const int start, const int stop) { mIndicesMother = {start, stop}; }
-  void setIndicesDaughter(const int start, const int stop) { mIndicesDaughter = {start, stop}; }
-  void setVxVyVz(const float vx, const float vy, const float vz)
-  {
-    mVx = vx;
-    mVy = vy;
-    mVz = vz;
-  }
-  void setPxPyPzE(const float px, const float py, const float pz, const float e)
-  {
-    mPx = px;
-    mPy = py;
-    mPz = pz;
-    mE = e;
-  }
-
-  // Getters
-  int pdgCode() const { return mPdgCode; }
-  int globalIndex() const { return mGlobalIndex; }
-  int collisionId() const { return mCollisionId; }
-  bool isAlive() const { return mIsAlive; }
-  bool isPrimary() const { return mIsPrimary; }
-  bool isFromMcParticles() const { return mIsFromMcParticles; }
-  float weight() const
-  {
-    static constexpr float Weight = 1.f;
-    return Weight;
-  }
-  uint8_t flags() const
-  {
-    static constexpr uint8_t Flags = 1;
-    return Flags;
-  }
-  float vt() const
-  {
-    static constexpr float Vt = 1.f;
-    return Vt;
-  }
-  int statusCode() const
-  {
-    static constexpr int StatusCode = 1;
-    return StatusCode;
-  }
-  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; }
-  float radius() const { return std::hypot(mVx, mVy); }
-  float r() const { return radius(); }
-  float pt() const { return std::hypot(mPx, mPy); }
-  float p() const { return std::hypot(mPx, mPy, mPz); }
-  float phi() const { return o2::constants::math::PI + std::atan2(-1.0f * py(), -1.0f * px()); }
-  float eta() const
-  {
-    // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943
-    static constexpr float Tolerance = 1e-7f;
-    if ((p() - mPz) < Tolerance) {
-      return (mPz < 0.0f) ? -100.0f : 100.0f;
-    } else {
-      return 0.5f * std::log((p() + mPz) / (p() - mPz));
-    }
-  }
-
-  float y() const
-  {
-    // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1922
-    static constexpr float Tolerance = 1e-7f;
-    if ((e() - mPz) < Tolerance) {
-      return (mPz < 0.0f) ? -100.0f : 100.0f;
-    } else {
-      return 0.5f * std::log((mE + mPz) / (mE - mPz));
-    }
-  }
-
-  bool hasDaughters() const { return (mIndicesDaughter[0] > 0); }
-  bool hasMothers() const { return (mIndicesMother[0] > 0); }
-  int getMotherIndexStart() const { return mIndicesMother[0]; }
-  int getMotherIndexStop() const { return mIndicesMother[1]; }
-  int getDaughterIndexStart() const { return mIndicesDaughter[0]; }
-  int getDaughterIndexStop() const { return mIndicesDaughter[1]; }
-  std::array<int, 2> getMothers() const { return mIndicesMother; }
-  std::array<int, 2> getDaughters() const { return mIndicesDaughter; }
-  std::span<const int> getMotherSpan() const { return hasMothers() ? std::span<const int>(mIndicesMother.data(), 2) : std::span<const int>(); }
-
-  bool hasNaN() const
-  {
-    return std::isnan(mPx) || std::isnan(mPy) || std::isnan(mPz) || std::isnan(mE) ||
-           std::isnan(mVx) || std::isnan(mVy) || std::isnan(mVz);
-  }
-
-  bool hasIndex() const
-  {
-    return (mGlobalIndex != -1);
-  }
-
- private:
-  int mPdgCode{}, mGlobalIndex{-1};
-  int mCollisionId{};
-  float mE{};
-  float mVx{}, mVy{}, mVz{};
-  float mPx{}, mPy{}, mPz{};
-  bool mIsAlive{}, mIsFromMcParticles{false};
-  bool mIsPrimary{};
-
-  std::array<int, 2> mIndicesMother{-1, -1}, mIndicesDaughter{-1, -1};
-};
-
 /// Function to convert a TLorentzVector into a perfect Track
 /// \param charge particle charge (integer)
 /// \param particle the particle to convert (TLorentzVector)
@@ -241,6 +104,17 @@ o2::track::TrackParCov convertMCParticleToO2Track(McParticleType& particle,
   return o2track;
 }
 
+/// returns velocity in centimeters per picoseconds
+/// \param momentum the momentum of the track
+/// \param mass the mass of the particle
+float computeParticleVelocity(float momentum, float mass);
+
+/// function to calculate track length of this track up to a certain radius
+/// \param track the input track (TrackParCov)
+/// \param radius the radius of the layer you're calculating the length to
+/// \param magneticField the magnetic field to use when propagating
+float computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField);
+
 } // namespace o2::upgrade
 
 #endif // ALICE3_CORE_TRACKUTILITIES_H_
diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx
index 1d1a301cd3f..0b69cba224b 100644
--- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx
+++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx
@@ -76,6 +76,7 @@ struct OnTheFlyDecayer {
   o2::upgrade::Decayer decayer;
   Service<o2::framework::O2DatabasePDG> pdgDB;
   std::map<int, std::vector<o2::upgrade::OTFParticle>> mDecayDaughters;
+  HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
 
   Configurable<int> seed{"seed", 0, "Set seed for particle decayer"};
   Configurable<float> magneticField{"magneticField", 20., "Magnetic field (kG)"};
@@ -83,9 +84,9 @@ struct OnTheFlyDecayer {
                                                 {DefaultParameters[0], NumDecays, NumParameters, ParticleNames, ParameterNames},
                                                 "Enable option for particle to be decayed: 0 - no, 1 - yes"};
 
-  HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
-
+  static constexpr float PicoToNano = 1.e-3f;
   int mCollisionId{-1};
+
   std::vector<int> mEnabledDecays;
   void init(o2::framework::InitContext&)
   {
@@ -132,12 +133,29 @@ struct OnTheFlyDecayer {
 
       particle.setIsAlive(false);
       std::vector<o2::upgrade::OTFParticle> decayStack = decayer.decayParticle(pdgDB, particle);
+      const float decayRadius = decayer.getDecayRadius();
+      const float trackVelocity = o2::upgrade::computeParticleVelocity(particle.p(), pdgDB->GetParticle(particle.pdgCode())->Mass());
+      const int charge = pdgDB->GetParticle(particle.pdgCode())->Charge() / 3;
+      float trackLength{-1.f};
+      if (!charge) {
+        const float dx = particle.vx() - decayer.getSecondaryVertexX();
+        const float dy = particle.vy() - decayer.getSecondaryVertexY();
+        const float dz = particle.vz() - decayer.getSecondaryVertexZ();
+        trackLength = std::hypot(dx, dy, dz);
+      } else {
+        o2::track::TrackParCov o2track;
+        o2::upgrade::convertOTFParticleToO2Track(particle, o2track, pdgDB);
+        trackLength = o2::upgrade::computeTrackLength(o2track, decayRadius, magneticField);
+      }
+
+      const float trackTimeNS = trackLength / trackVelocity * PicoToNano;
       particle.setIndicesDaughter(allParticles.size(), allParticles.size() + (decayStack.size() - 1));
       for (o2::upgrade::OTFParticle daughter : decayStack) {
         daughter.setIndicesMother(i, i);
         daughter.setCollisionId(mCollisionId);
         daughter.setIsAlive(true);
         daughter.setIsPrimary(false);
+        daughter.setProductionTime(trackTimeNS);
         allParticles.push_back(daughter);
         ndau++;
       }
@@ -174,7 +192,7 @@ struct OnTheFlyDecayer {
         histos.fill(HIST("hNaNBookkeeping"), 0);
       }
 
-      // todo: status codes and vt
+      // todo: status codes
       tableMcParticlesWithDau(otfParticle.collisionId(), otfParticle.pdgCode(), otfParticle.statusCode(),
                               otfParticle.flags(), otfParticle.getMotherSpan(), otfParticle.getDaughters().data(), otfParticle.weight(),
                               otfParticle.px(), otfParticle.py(), otfParticle.pz(), otfParticle.e(),
diff --git a/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx
index b327b098a06..4266054c81d 100644
--- a/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx
+++ b/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx
@@ -448,81 +448,6 @@ struct OnTheFlyTofPid {
     return outerTOFLayer.isTrackInActiveArea(track);
   }
 
-  /// function to calculate track length of this track up to a certain radius
-  /// \param track the input track
-  /// \param radius the radius of the layer you're calculating the length to
-  /// \param magneticField the magnetic field to use when propagating
-  static float computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField)
-  {
-    // don't make use of the track parametrization
-    float length = -100;
-
-    o2::math_utils::CircleXYf_t trcCircle;
-    float sna, csa;
-    track.getCircleParams(magneticField, trcCircle, sna, csa);
-
-    // distance between circle centers (one circle is at origin -> easy)
-    const float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC);
-
-    // condition of circles touching - if not satisfied returned length will be -100
-    if (centerDistance < trcCircle.rC + radius && centerDistance > std::fabs(trcCircle.rC - radius)) {
-      length = 0.0f;
-
-      // base radical direction
-      const float ux = trcCircle.xC / centerDistance;
-      const float uy = trcCircle.yC / centerDistance;
-      // calculate perpendicular vector (normalized) for +/- displacement
-      const float vx = -uy;
-      const float vy = +ux;
-      // calculate coordinate for radical line
-      const float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + radius * radius) / (2.0f * centerDistance);
-      // calculate absolute displacement from center-to-center axis
-      const float displace = (0.5f / centerDistance) * std::sqrt(
-                                                         (-centerDistance + trcCircle.rC - radius) *
-                                                         (-centerDistance - trcCircle.rC + radius) *
-                                                         (-centerDistance + trcCircle.rC + radius) *
-                                                         (centerDistance + trcCircle.rC + radius));
-
-      // possible intercept points of track and TOF layer in 2D plane
-      const float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy};
-      const float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy};
-
-      // decide on correct intercept point
-      std::array<float, 3> mom;
-      track.getPxPyPzGlo(mom);
-      const float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1];
-      const float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1];
-
-      // get start point
-      std::array<float, 3> startPoint;
-      track.getXYZGlo(startPoint);
-
-      float cosAngle = -1000, modulus = -1000;
-
-      if (scalarProduct1 > scalarProduct2) {
-        modulus = std::hypot(point1[0] - trcCircle.xC, point1[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
-        cosAngle = (point1[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point1[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC);
-      } else {
-        modulus = std::hypot(point2[0] - trcCircle.xC, point2[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
-        cosAngle = (point2[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point2[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC);
-      }
-      cosAngle /= modulus;
-      length = trcCircle.rC * std::acos(cosAngle);
-      length *= std::sqrt(1.0f + track.getTgl() * track.getTgl());
-    }
-    return length;
-  }
-
-  /// returns velocity in centimeters per picoseconds
-  /// \param momentum the momentum of the tarck
-  /// \param mass the mass of the particle
-  float computeParticleVelocity(float momentum, float mass)
-  {
-    const float a = momentum / mass;
-    // uses light speed in cm/ps so output is in those units
-    return o2::constants::physics::LightSpeedCm2PS * a / std::sqrt((1.f + a * a));
-  }
-
   struct TracksWithTime {
     TracksWithTime(int pdgCode,
                    std::pair<float, float> innerTOFTime,
@@ -696,8 +621,8 @@ struct OnTheFlyTofPid {
       }
       float trackLengthInnerTOF = -1, trackLengthOuterTOF = -1;
       if (xPv > kTrkXThreshold) {
-        trackLengthInnerTOF = computeTrackLength(o2track, simConfig.innerTOFRadius, mMagneticField);
-        trackLengthOuterTOF = computeTrackLength(o2track, simConfig.outerTOFRadius, mMagneticField);
+        trackLengthInnerTOF = o2::upgrade::computeTrackLength(o2track, simConfig.innerTOFRadius, mMagneticField);
+        trackLengthOuterTOF = o2::upgrade::computeTrackLength(o2track, simConfig.outerTOFRadius, mMagneticField);
       }
 
       // Check if the track hit a sensitive area of the TOF
@@ -726,7 +651,7 @@ struct OnTheFlyTofPid {
         upgradeTofMC(-999.f, -999.f, -999.f, -999.f);
         continue;
       }
-      const float v = computeParticleVelocity(o2track.getP(), pdgInfo->Mass());
+      const float v = o2::upgrade::computeParticleVelocity(o2track.getP(), pdgInfo->Mass());
       const float expectedTimeInnerTOF = trackLengthInnerTOF > 0 ? trackLengthInnerTOF / v + eventCollisionTimePS : -999.f; // arrival time to the Inner TOF in ps
       const float expectedTimeOuterTOF = trackLengthOuterTOF > 0 ? trackLengthOuterTOF / v + eventCollisionTimePS : -999.f; // arrival time to the Outer TOF in ps
       upgradeTofMC(expectedTimeInnerTOF, trackLengthInnerTOF, expectedTimeOuterTOF, trackLengthOuterTOF);
@@ -744,8 +669,8 @@ struct OnTheFlyTofPid {
         xPv = recoTrack.getX();
       }
       if (xPv > kTrkXThreshold) {
-        trackLengthRecoInnerTOF = computeTrackLength(recoTrack, simConfig.innerTOFRadius, mMagneticField);
-        trackLengthRecoOuterTOF = computeTrackLength(recoTrack, simConfig.outerTOFRadius, mMagneticField);
+        trackLengthRecoInnerTOF = o2::upgrade::computeTrackLength(recoTrack, simConfig.innerTOFRadius, mMagneticField);
+        trackLengthRecoOuterTOF = o2::upgrade::computeTrackLength(recoTrack, simConfig.outerTOFRadius, mMagneticField);
       }
 
       // cache the track info needed for the event time calculation
@@ -866,7 +791,7 @@ struct OnTheFlyTofPid {
         nSigmaOuterTOF[ii] = -100;
 
         momentumHypotheses[ii] = rigidity * kParticleCharges[ii]; // Total momentum for this hypothesis
-        const float v = computeParticleVelocity(momentumHypotheses[ii], kParticleMasses[ii]);
+        const float v = o2::upgrade::computeParticleVelocity(momentumHypotheses[ii], kParticleMasses[ii]);
 
         expectedTimeInnerTOF[ii] = trackLengthInnerTOF / v;
         expectedTimeOuterTOF[ii] = trackLengthOuterTOF / v;