[ALICE3] Add particle production time to otf decayer

ALICE3/Core/Decayer.h

@@ -28,7 +28,7 @@
 
 #include <TDecayChannel.h> // IWYU pragma: keep
 #include <TGenPhaseSpace.h>
 #include <TLorentzVector.h>
 #include <TRandom3.h>
 
 #include <cmath>
@@ -61,28 +61,27 @@
     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 {
       o2::track::TrackParCov track;
       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);
@@ -115,7 +114,7 @@
       return {};
     }
 
     TLorentzVector tlv(px, py, particle.pz(), e);
     TGenPhaseSpace decay;
     decay.SetDecay(tlv, dauMasses.size(), dauMasses.data());
     decay.Generate();
@@ -123,27 +122,34 @@
     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.setVxVyVz(mVx, mVy, mVz);
       particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E());
       decayProducts.push_back(particle);
     }
 
     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

ALICE3/Core/TrackUtilities.cxx

@@ -20,14 +20,14 @@
 #include <MathUtils/Utils.h>
 #include <ReconstructionDataFormats/Track.h>
 
 #include <TLorentzVector.h>
 
 #include <array>
 #include <cmath>
 #include <vector>
 
 void o2::upgrade::convertTLorentzVectorToO2Track(const int charge,
                                                  const TLorentzVector particle,
                                                  const std::vector<double> productionVertex,
                                                  o2::track::TrackParCov& o2track)
 {
@@ -45,3 +45,71 @@
   // 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;
+}

ALICE3/Core/TrackUtilities.h

@@ -22,7 +22,7 @@
 
 #include <ReconstructionDataFormats/Track.h>
 
 #include <TLorentzVector.h>
 
 #include <array>
 #include <cmath>
@@ -37,8 +37,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);
 
@@ -49,7 +49,7 @@
 /// \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,
@@ -104,6 +104,17 @@
   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_

ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx

@@ -77,16 +77,17 @@ 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)"};
   Configurable<LabeledArray<int>> enabledDecays{"enabledDecays",
                                                 {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&)
   {
@@ -133,12 +134,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++;
       }
@@ -175,7 +193,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(),

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
@@ -730,7 +655,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);
@@ -748,8 +673,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
@@ -870,7 +795,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;