Merge branch 'master' into const-PWGCF

Author: vkucera

ALICE3/Core/Decayer.h

@@ -31,13 +31,12 @@
 #include <TLorentzVector.h>
 #include <TRandom3.h>
 
+#include <array>
 #include <cmath>
 #include <cstddef>
 #include <vector>
 
-namespace o2
-{
-namespace upgrade
+namespace o2::upgrade
 {
 
 class Decayer
@@ -47,44 +46,27 @@ class Decayer
   Decayer() = default;
 
   template <typename TDatabase>
-  std::vector<o2::upgrade::OTFParticle> decayParticle(const TDatabase& pdgDB, const OTFParticle& particle)
+  std::vector<o2::upgrade::OTFParticle> decayParticle(const OTFParticle& particle, const TDatabase& pdgDB)
   {
-    const auto& particleInfo = pdgDB->GetParticle(particle.pdgCode());
+    auto particleInfo = pdgDB->GetParticle(particle.pdgCode());
     if (!particleInfo) {
       return {};
     }
 
     const int charge = particleInfo->Charge() / 3;
     const double mass = particleInfo->Mass();
-
-    const double u = mRand3.Uniform(0.001, 0.999);
-    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 px, py, e;
+    std::array<double, 3> decayVtx = generateDecayVertex<double>(particle, pdgDB);
+    mVx = decayVtx[0];
+    mVy = decayVtx[1];
+    mVz = decayVtx[2];
+    double px{}, py{}, e{};
 
     if (!charge) {
-      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{};
-      track.getCircleParams(mBz, circle, sna, csa);
-      const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl());
-      const double theta = rxy / circle.rC;
-
-      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);
+      px = particle.px() * std::cos(mTheta) - particle.py() * std::sin(mTheta);
+      py = particle.py() * std::cos(mTheta) + particle.px() * std::sin(mTheta);
     }
 
     double brTotal = 0.;
@@ -133,6 +115,42 @@ class Decayer
     return decayProducts;
   }
 
+  template <typename T = float, typename TDatabase, typename TParticle>
+  std::array<T, 3> generateDecayVertex(const TParticle& particle, const TDatabase& pdgDB)
+  {
+    std::array<T, 3> decayVertex{};
+    auto particleInfo = pdgDB->GetParticle(particle.pdgCode());
+    if (!particleInfo) {
+      return {};
+    }
+
+    const int charge = particleInfo->Charge() / 3;
+    const double mass = particleInfo->Mass();
+    const double u = mRand3.Uniform(0.001, 0.999);
+    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);
+
+    if (!charge) {
+      decayVertex[0] = particle.vx() + rxyz * (particle.px() / particle.p());
+      decayVertex[1] = particle.vy() + rxyz * (particle.py() / particle.p());
+      decayVertex[2] = particle.vz() + rxyz * (particle.pz() / particle.p());
+    } else {
+      o2::math_utils::CircleXYf_t circle;
+      o2::track::TrackParCov track = o2::upgrade::convertMCParticleToO2Track(particle, pdgDB);
+
+      float sna{}, csa{};
+      track.getCircleParams(mBz, circle, sna, csa);
+      const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl());
+      mTheta = rxy / circle.rC;
+
+      decayVertex[0] = ((particle.vx() - circle.xC) * std::cos(mTheta) - (particle.vy() - circle.yC) * std::sin(mTheta)) + circle.xC;
+      decayVertex[1] = ((particle.vy() - circle.yC) * std::cos(mTheta) + (particle.vx() - circle.xC) * std::sin(mTheta)) + circle.yC;
+      decayVertex[2] = particle.vz() + rxyz * (particle.pz() / track.getP());
+    }
+    return decayVertex;
+  }
+
   // Setters
   void setBField(const double b) { mBz = b; }
   void setSeed(const int seed)
@@ -142,18 +160,18 @@ class Decayer
   }
 
   // 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)); }
+  [[nodiscard]] float getSecondaryVertexX() const { return static_cast<float>(mVx); }
+  [[nodiscard]] float getSecondaryVertexY() const { return static_cast<float>(mVy); }
+  [[nodiscard]] float getSecondaryVertexZ() const { return static_cast<float>(mVz); }
+  [[nodiscard]] float getDecayRadius() const { return static_cast<float>(std::hypot(mVx, mVy)); }
 
  private:
   double mBz{20.}; // kG
   double mVx{-1.}, mVy{-1.}, mVz{-1.};
-  TRandom3 mRand3{};
+  double mTheta{};
+  TRandom3 mRand3;
 };
 
-} // namespace upgrade
-} // namespace o2
+} // namespace o2::upgrade
 
 #endif // ALICE3_CORE_DECAYER_H_

ALICE3/Core/DelphesO2LutWriter.cxx

@@ -376,7 +376,7 @@ void DelphesO2LutWriter::diagonalise(lutEntry_t& lutEntry)
   // m.Print();
   TMatrixDSymEigen eigen(m);
   // eigenvalues vector
-  TVectorD eigenVal = eigen.GetEigenValues();
+  const TVectorD& eigenVal = eigen.GetEigenValues();
   for (int i = 0; i < kEig; ++i)
     lutEntry.eigval[i] = eigenVal[i];
   // eigenvectors matrix

ALICE3/Core/FastTracker.cxx

@@ -51,7 +51,7 @@ namespace fastsim
 
 // +-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+-~-<*>-~-+
 
-DetLayer* FastTracker::AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi, float resZ, float eff, int type)
+DetLayer* FastTracker::AddLayer(const TString& name, float r, float z, float x0, float xrho, float resRPhi, float resZ, float eff, int type)
 {
   LOG(debug) << "Adding layer " << name << " r=" << r << " z=" << z << " x0=" << x0 << " xrho=" << xrho << " resRPhi=" << resRPhi << " resZ=" << resZ << " eff=" << eff << " type=" << type;
   DetLayer newLayer(name, r, z, x0, xrho, resRPhi, resZ, eff, type);
@@ -154,7 +154,7 @@ void FastTracker::AddTPC(float phiResMean, float zResMean)
   }
 }
 
-void FastTracker::AddGenericDetector(o2::fastsim::GeometryEntry configMap, o2::ccdb::BasicCCDBManager* ccdbManager)
+void FastTracker::AddGenericDetector(const o2::fastsim::GeometryEntry& configMap, o2::ccdb::BasicCCDBManager* ccdbManager)
 {
   // Layers
   for (const auto& layer : configMap.getLayerNames()) {
@@ -586,7 +586,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa
   m.SetMatrixArray(reinterpret_cast<double*>(fcovm));
   TMatrixDSymEigen eigen(m);
   TMatrixD eigVec = eigen.GetEigenVectors();
-  TVectorD eigVal = eigen.GetEigenValues();
+  const TVectorD& eigVal = eigen.GetEigenValues();
   bool negEigVal = false;
   for (int ii = 0; ii < 5; ii++) {
     if (eigVal[ii] < 0.0f)

ALICE3/Core/FastTracker.h

@@ -46,7 +46,7 @@ class FastTracker
   virtual ~FastTracker() {}
 
   // Layer and layer configuration
-  DetLayer* AddLayer(TString name, float r, float z, float x0, float xrho, float resRPhi = 0.0f, float resZ = 0.0f, float eff = 0.0f, int type = 0);
+  DetLayer* AddLayer(const TString& name, float r, float z, float x0, float xrho, float resRPhi = 0.0f, float resZ = 0.0f, float eff = 0.0f, int type = 0);
 
   /// Add a dead region in phi for a specific layer
   /// \param layerName Name of the layer to modify
@@ -59,11 +59,11 @@ class FastTracker
   size_t GetNLayers() const { return layers.size(); }
   bool IsLayerInert(const int layer) const { return layers[layer].isInert(); }
   void ClearLayers() { layers.clear(); }
-  void SetRadiationLength(const std::string layerName, float x0) { layers[GetLayerIndex(layerName)].setRadiationLength(x0); }
-  void SetRadius(const std::string layerName, float r) { layers[GetLayerIndex(layerName)].setRadius(r); }
-  void SetResolutionRPhi(const std::string layerName, float resRPhi) { layers[GetLayerIndex(layerName)].setResolutionRPhi(resRPhi); }
-  void SetResolutionZ(const std::string layerName, float resZ) { layers[GetLayerIndex(layerName)].setResolutionZ(resZ); }
-  void SetResolution(const std::string layerName, float resRPhi, float resZ)
+  void SetRadiationLength(const std::string& layerName, float x0) { layers[GetLayerIndex(layerName)].setRadiationLength(x0); }
+  void SetRadius(const std::string& layerName, float r) { layers[GetLayerIndex(layerName)].setRadius(r); }
+  void SetResolutionRPhi(const std::string& layerName, float resRPhi) { layers[GetLayerIndex(layerName)].setResolutionRPhi(resRPhi); }
+  void SetResolutionZ(const std::string& layerName, float resZ) { layers[GetLayerIndex(layerName)].setResolutionZ(resZ); }
+  void SetResolution(const std::string& layerName, float resRPhi, float resZ)
   {
     SetResolutionRPhi(layerName, resRPhi);
     SetResolutionZ(layerName, resZ);
@@ -80,7 +80,7 @@ class FastTracker
    *
    * @param configMap Configuration map describing the detector.
    */
-  void AddGenericDetector(o2::fastsim::GeometryEntry configMap, o2::ccdb::BasicCCDBManager* ccdbManager = nullptr);
+  void AddGenericDetector(const o2::fastsim::GeometryEntry& configMap, o2::ccdb::BasicCCDBManager* ccdbManager = nullptr);
 
   void Print();
 

ALICE3/Core/FlatLutWriter.cxx

@@ -392,7 +392,7 @@ void FlatLutWriter::diagonalise(lutEntry_t& lutEntry)
   TMatrixDSymEigen eigen(m);
 
   // Eigenvalues
-  TVectorD eigenVal = eigen.GetEigenValues();
+  const TVectorD& eigenVal = eigen.GetEigenValues();
   for (int i = 0; i < kEig; ++i)
     lutEntry.eigval[i] = eigenVal[i];
 

ALICE3/Core/GeometryContainer.cxx

@@ -159,7 +159,7 @@ void GeometryEntry::replaceValue(const std::string& layerName, const std::string
   setValue(layerName, key, value);
 }
 
-std::string GeometryEntry::accessFile(const std::string& path, const std::string downloadPath, o2::ccdb::BasicCCDBManager* ccdb, int timeoutSeconds)
+std::string GeometryEntry::accessFile(const std::string& path, const std::string& downloadPath, o2::ccdb::BasicCCDBManager* ccdb, int timeoutSeconds)
 {
 
   if (path.rfind("ccdb:", 0) == 0) {

ALICE3/Core/GeometryContainer.h

@@ -33,7 +33,7 @@ namespace o2::fastsim
 struct GeometryEntry {
   // Default constructor
   GeometryEntry() = default;
-  explicit GeometryEntry(std::string filename, o2::ccdb::BasicCCDBManager* ccdb = nullptr)
+  explicit GeometryEntry(const std::string& filename, o2::ccdb::BasicCCDBManager* ccdb = nullptr)
   {
     mFileName = accessFile(filename, "./.ALICE3/Configuration/", ccdb);
     mConfigurations = GeometryEntry::parseTEnvConfiguration(mFileName, mLayerNames);
@@ -59,7 +59,7 @@ struct GeometryEntry {
    * @param timeoutSeconds If positive, then this function will wait for these seconds after download before removing the downloaded file.
    * @return The local path to the file, either the original local path or the path to the retrieved file from ccdb
    */
-  static std::string accessFile(const std::string& path, const std::string downloadPath = "/tmp/GeometryContainer/", o2::ccdb::BasicCCDBManager* ccdb = nullptr, int timeoutSeconds = 0);
+  static std::string accessFile(const std::string& path, const std::string& downloadPath = "/tmp/GeometryContainer/", o2::ccdb::BasicCCDBManager* ccdb = nullptr, int timeoutSeconds = 0);
 
   std::map<std::string, std::map<std::string, std::string>> getConfigurations() const { return mConfigurations; }
   std::map<std::string, std::string> getConfiguration(const std::string& layerName) const;

ALICE3/DataModel/OTFTOF.h

@@ -150,6 +150,26 @@ DECLARE_SOA_DYNAMIC_COLUMN(NSigmaOuterTOF, nSigmaOuterTOF, //! General function
                              }
                            });
 
+DECLARE_SOA_COLUMN(NSigmaSplusInnerTOF, nSigmaSplusInnerTOF, float);   //! NSigma Sigma plus InnerTOF
+DECLARE_SOA_COLUMN(NSigmaSminusInnerTOF, nSigmaSminusInnerTOF, float); //! NSigma Sigma minus InnerTOF
+DECLARE_SOA_COLUMN(NSigmaXiInnerTOF, nSigmaXiInnerTOF, float);         //! NSigma Xi InnerTOF
+DECLARE_SOA_COLUMN(NSigmaOmegaInnerTOF, nSigmaOmegaInnerTOF, float);   //! NSigma Omega InnerTOF
+
+DECLARE_SOA_COLUMN(NSigmaSplusOuterTOF, nSigmaSplusOuterTOF, float);   //! NSigma Sigma plus OuterTOF
+DECLARE_SOA_COLUMN(NSigmaSminusOuterTOF, nSigmaSminusOuterTOF, float); //! NSigma Sigma minus OuterTOF
+DECLARE_SOA_COLUMN(NSigmaXiOuterTOF, nSigmaXiOuterTOF, float);         //! NSigma Xi OuterTOF
+DECLARE_SOA_COLUMN(NSigmaOmegaOuterTOF, nSigmaOmegaOuterTOF, float);   //! NSigma Omega OuterTOF
+
+DECLARE_SOA_COLUMN(InnerTOFExpectedTimeSp, innerTOFExpectedTimeSp, float); //! Reconstructed expected time at the InnerTOF for the Sigma plus mass hypotheses
+DECLARE_SOA_COLUMN(InnerTOFExpectedTimeSm, innerTOFExpectedTimeSm, float); //! Reconstructed expected time at the InnerTOF for the Sigma minus mass hypotheses
+DECLARE_SOA_COLUMN(InnerTOFExpectedTimeXi, innerTOFExpectedTimeXi, float); //! Reconstructed expected time at the InnerTOF for the Xi mass hypotheses
+DECLARE_SOA_COLUMN(InnerTOFExpectedTimeOm, innerTOFExpectedTimeOm, float); //! Reconstructed expected time at the InnerTOF for the Omega mass hypotheses
+
+DECLARE_SOA_COLUMN(OuterTOFExpectedTimeSp, outerTOFExpectedTimeSp, float); //! Reconstructed expected time at the OuterTOF for the Sigma plus mass hypotheses
+DECLARE_SOA_COLUMN(OuterTOFExpectedTimeSm, outerTOFExpectedTimeSm, float); //! Reconstructed expected time at the OuterTOF for the Sigma minus mass hypotheses
+DECLARE_SOA_COLUMN(OuterTOFExpectedTimeXi, outerTOFExpectedTimeXi, float); //! Reconstructed expected time at the OuterTOF for the Xi mass hypotheses
+DECLARE_SOA_COLUMN(OuterTOFExpectedTimeOm, outerTOFExpectedTimeOm, float); //! Reconstructed expected time at the OuterTOF for the Omega mass hypotheses
+
 } // namespace upgrade_tof
 
 DECLARE_SOA_TABLE(UpgradeTofMCs, "AOD", "UPGRADETOFMC",
@@ -222,9 +242,28 @@ DECLARE_SOA_TABLE(UpgradeTofExpectedTimes, "AOD", "UPGRADETOFEXPT",
                   upgrade_tof::OuterTOFExpectedTimeHe3,
                   upgrade_tof::OuterTOFExpectedTimeAl);
 
+DECLARE_SOA_TABLE(UpgradeTofShortLiveds, "AOD", "UPGRTOFSHRTLVD",
+                  upgrade_tof::NSigmaSplusInnerTOF,
+                  upgrade_tof::NSigmaSminusInnerTOF,
+                  upgrade_tof::NSigmaXiInnerTOF,
+                  upgrade_tof::NSigmaOmegaInnerTOF,
+                  upgrade_tof::NSigmaSplusOuterTOF,
+                  upgrade_tof::NSigmaSminusOuterTOF,
+                  upgrade_tof::NSigmaXiOuterTOF,
+                  upgrade_tof::NSigmaOmegaOuterTOF,
+                  upgrade_tof::InnerTOFExpectedTimeSp,
+                  upgrade_tof::InnerTOFExpectedTimeSm,
+                  upgrade_tof::InnerTOFExpectedTimeXi,
+                  upgrade_tof::InnerTOFExpectedTimeOm,
+                  upgrade_tof::OuterTOFExpectedTimeSp,
+                  upgrade_tof::OuterTOFExpectedTimeSm,
+                  upgrade_tof::OuterTOFExpectedTimeXi,
+                  upgrade_tof::OuterTOFExpectedTimeOm);
+
 using UpgradeTofMC = UpgradeTofMCs::iterator;
 using UpgradeTof = UpgradeTofs::iterator;
 using UpgradeTofExpectedTime = UpgradeTofExpectedTimes::iterator;
+using UpgradeTofShortLived = UpgradeTofShortLiveds::iterator;
 
 } // namespace o2::aod