[ALICE3] add production time to daughters in otf-decayer

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>
@@ -27,7 +28,7 @@
 
 #include <TDecayChannel.h> // IWYU pragma: keep
 #include <TGenPhaseSpace.h>
 #include <TLorentzVector.h>
 #include <TRandom3.h>
 
 #include <cmath>
@@ -60,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);
@@ -114,7 +114,7 @@
       return {};
     }
 
     TLorentzVector tlv(px, py, particle.pz(), e);
     TGenPhaseSpace decay;
     decay.SetDecay(tlv, dauMasses.size(), dauMasses.data());
     decay.Generate();
@@ -122,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/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_

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;
+}