Merge branch 'dev' into dev

Author: njacazio

DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackFwd.h

@@ -31,6 +31,12 @@ using SMatrix55Sym = ROOT::Math::SMatrix<double, 5, 5, ROOT::Math::MatRepSym<dou
 using SMatrix55Std = ROOT::Math::SMatrix<double, 5>;
 using SMatrix5 = ROOT::Math::SVector<Double_t, 5>;
 
+template <typename value_T>
+class TrackParametrization; // fwd declaration for conversion method
+
+template <typename value_T>
+class TrackParametrizationWithError; // fwd declaration for conversion method
+
 class TrackParFwd
 { // Forward track parameterization, kinematics only.
  public:
@@ -42,6 +48,9 @@ class TrackParFwd
   TrackParFwd(TrackParFwd&&) = delete;
   TrackParFwd& operator=(TrackParFwd&&) = delete;
 
+  template <typename T>
+  void toBarrelTrackPar(TrackParametrization<T>& t) const;
+
   /// return Z coordinate (cm)
   Double_t getZ() const { return mZ; }
   /// set Z coordinate (cm)
@@ -145,6 +154,9 @@ class TrackParCovFwd : public TrackParFwd
   TrackParCovFwd& operator=(const TrackParCovFwd& tpf) = default;
   TrackParCovFwd(const Double_t z, const SMatrix5& parameters, const SMatrix55Sym& covariances, const Double_t chi2);
 
+  template <typename T>
+  void toBarrelTrackParCov(TrackParametrizationWithError<T>& t) const;
+
   const SMatrix55Sym& getCovariances() const { return mCovariances; }
   void setCovariances(const SMatrix55Sym& covariances) { mCovariances = covariances; }
   void deleteCovariances() { mCovariances = SMatrix55Sym(); }

DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h

@@ -66,6 +66,9 @@ class DCA;
 
 namespace track
 {
+
+class TrackParFwd; // fwd declaration for conversion method
+
 // aliases for track elements
 enum ParLabels : int { kY,
                        kZ,
@@ -252,6 +255,7 @@ class TrackParametrization
   GPUd() void printParam() const;
   GPUd() void printParamHexadecimal();
 #ifndef GPUCA_ALIGPUCODE
+  void toFwdTrackPar(TrackParFwd& t) const;
   std::string asString() const;
   std::string asStringHexadecimal();
   size_t hash() const { return hash(getX(), getAlpha(), getY(), getZ(), getSnp(), getTgl(), getQ2Pt()); }

DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h

@@ -25,6 +25,8 @@ namespace o2
 namespace track
 {
 
+class TrackParCovFwd; // fwd declaration for conversion method
+
 template <typename value_T = float>
 class TrackParametrizationWithError : public TrackParametrization<value_T>
 { // track+error parameterization
@@ -82,6 +84,7 @@ class TrackParametrizationWithError : public TrackParametrization<value_T>
   GPUd() void print() const;
   GPUd() void printHexadecimal();
 #ifndef GPUCA_ALIGPUCODE
+  bool toFwdTrackParCov(TrackParCovFwd& t) const;
   std::string asString() const;
   std::string asStringHexadecimal();
 #endif

DataFormats/Reconstruction/src/TrackFwd.cxx

@@ -10,6 +10,8 @@
 // or submit itself to any jurisdiction.
 
 #include "ReconstructionDataFormats/TrackFwd.h"
+#include "ReconstructionDataFormats/TrackParametrization.h"
+#include "ReconstructionDataFormats/TrackParametrizationWithError.h"
 #include "Math/MatrixFunctions.h"
 #include <GPUCommonLogger.h>
 
@@ -41,7 +43,7 @@ void TrackParFwd::propagateParamToZlinear(double zEnd)
   auto dZ = (zEnd - getZ());
   auto phi0 = getPhi();
   auto [sinphi0, cosphi0] = o2::math_utils::sincosd(phi0);
-  auto invtanl0 = 1.0 / getTanl();
+  auto invtanl0 = 1.0 / getTgl();
   auto n = dZ * invtanl0;
   mParameters(0) += n * cosphi0;
   mParameters(1) += n * sinphi0;
@@ -572,5 +574,79 @@ void TrackParCovFwd::propagateToDCAhelix(double zField, const std::array<double,
   return;
 }
 
+template <typename T>
+void TrackParFwd::toBarrelTrackPar(TrackParametrization<T>& t) const
+{
+  // we select the barrel frame with alpha = phi, then by construction the snp is 0
+  auto alpha = getPhi();
+  auto csa = TMath::Cos(alpha), sna = TMath::Sin(alpha);
+  t.setAlpha(alpha);
+  t.setX(csa * getX() + sna * getY());
+  t.setY(-sna * getX() + csa * getY());
+  t.setZ(getZ());
+  t.setSnp(0);
+  t.setTgl(getTanl());
+  t.setQ2Pt(getInvQPt());
+}
+
+template <typename T>
+void TrackParCovFwd::toBarrelTrackParCov(TrackParametrizationWithError<T>& t) const
+{
+  // We select the barrel frame with alpha = phi, then by construction the snp is 0
+  auto alpha = getPhi();
+  auto csa = TMath::Cos(alpha), sna = TMath::Sin(alpha);
+  t.setAlpha(alpha);
+  t.setX(csa * getX() + sna * getY());
+  t.setY(-sna * getX() + csa * getY());
+  t.setZ(getZ());
+  t.setSnp(0);
+  t.setTgl(getTgl());
+  t.setQ2Pt(getInvQPt());
+  /*
+    The standard Jacobian d{barrel_param} / d{fwd_param} should be augmented by the "forward" uncertainty
+    in X_barrel translated to uncertainty in Z, i.e:
+    A fwd param variation delta_x, delta_y leads to barrel frame coordinate variaion
+    delta_xb =  csa delta_x + sna delta_y
+    delta_yb = -sna delta_x + csa delta_y
+    with dx_b/dz = csp/tgL = 1/tgL, dy_b/dz = snp/tgL = 0  (for phi 0 in the tracking frame) the variation of delta_xb would require
+    a shift in delta_zb = -tgL delta_xb to stat at the same X.
+    So, for alpha=phi (-> snp=0) choice the full Jacobian fwd->barrel is:
+    / -sna       csa      0 0 0  \
+    | -tgL*csa  -tgL*sna  0 0 0  |
+    |  0          0       1 0 0  |
+    |  0          0       0 1 0  |
+    \  0          0       0 0 1 /
+  */
+  const T a1 = -sna;
+  const T a2 = csa;
+  const T b1 = -getTgl() * csa;
+  const T b2 = -getTgl() * sna;
+  const T cphi = 1;
+  const auto& C = getCovariances();
+  typename TrackParametrizationWithError<T>::covMat_t covBarrel = {
+    T(a1 * a1 * C(0, 0) + 2 * a1 * a2 * C(0, 1) + a2 * a2 * C(1, 1)),         // kSigY2
+    T(a1 * b1 * C(0, 0) + (a1 * b2 + a2 * b1) * C(0, 1) + a2 * b2 * C(1, 1)), // kSigZY
+    T(b1 * b1 * C(0, 0) + 2 * b1 * b2 * C(0, 1) + b2 * b2 * C(1, 1)),         // kSigZ2
+    T(csa * (a1 * C(0, 2) + a2 * C(1, 2))),                                   // kSigSnpY
+    T(csa * (b1 * C(0, 2) + b2 * C(1, 2))),                                   // kSigSnpZ
+    T(csa * csa * C(2, 2)),                                                   // kSigSnp2
+    T(a1 * C(0, 3) + a2 * C(1, 3)),                                           // kSigTglY
+    T(b1 * C(0, 3) + b2 * C(1, 3)),                                           // kSigTglZ
+    T(csa * C(2, 3)),                                                         // kSigTglSnp
+    T(C(3, 3)),                                                               // kSigTgl2
+    T(a1 * C(0, 4) + a2 * C(1, 4)),                                           // kSigQ2PtY
+    T(b1 * C(0, 4) + b2 * C(1, 4)),                                           // kSigQ2PtZ
+    T(csa * C(2, 4)),                                                         // kSigQ2PtSnp
+    T(C(3, 4)),                                                               // kSigQ2PtTgl
+    T(C(4, 4))                                                                // kSigQ2Pt2
+  };
+  t.setCov(covBarrel);
+}
+
+template void TrackParFwd::toBarrelTrackPar<float>(TrackParametrization<float>&) const;
+template void TrackParFwd::toBarrelTrackPar<double>(TrackParametrization<double>&) const;
+template void TrackParCovFwd::toBarrelTrackParCov<float>(TrackParametrizationWithError<float>&) const;
+template void TrackParCovFwd::toBarrelTrackParCov<double>(TrackParametrizationWithError<double>&) const;
+
 } // namespace track
 } // namespace o2

DataFormats/Reconstruction/src/TrackParametrization.cxx

@@ -26,6 +26,7 @@
 
 #ifndef GPUCA_ALIGPUCODE
 #include <fmt/printf.h>
+#include "ReconstructionDataFormats/TrackFwd.h"
 #endif
 
 using namespace o2::gpu;
@@ -985,6 +986,21 @@ GPUd() typename TrackParametrization<value_T>::value_t TrackParametrization<valu
   return ttmp.propagateParamToDCA({xmv, ymv, zmv}, b, &dca) ? dca[1] : -9999.;
 }
 
+#ifndef GPUCA_ALIGPUCODE
+//______________________________________________
+template <typename value_T>
+void TrackParametrization<value_T>::toFwdTrackPar(TrackParFwd& t) const
+{
+  auto p = getXYZGlo();
+  t.setZ(p.Z());
+  t.setX(p.X());
+  t.setY(p.Y());
+  t.setPhi(getPhi());
+  t.setTanl(getTgl());
+  t.setInvQPt(getQ2Pt());
+}
+#endif
+
 namespace o2::track
 {
 #if !defined(GPUCA_GPUCODE) || defined(GPUCA_GPUCODE_DEVICE) // FIXME: DR: WORKAROUND to avoid CUDA bug creating host symbols for device code.

DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx

@@ -21,6 +21,7 @@
 
 #ifndef GPUCA_ALIGPUCODE
 #include <fmt/printf.h>
+#include "ReconstructionDataFormats/TrackFwd.h"
 #endif
 
 using namespace o2::track;
@@ -1768,6 +1769,64 @@ GPUd() void TrackParametrizationWithError<value_T>::printHexadecimal()
 #endif
 }
 
+#ifndef GPUCA_ALIGPUCODE
+//______________________________________________________________
+template <typename value_T>
+bool TrackParametrizationWithError<value_T>::toFwdTrackParCov(TrackParCovFwd& t) const
+{
+  auto p = this->getXYZGlo();
+  t.setZ(p.Z());
+  t.setX(p.X());
+  t.setY(p.Y());
+  t.setPhi(this->getPhi());
+  t.setTanl(this->getTgl());
+  t.setInvQPt(this->getQ2Pt());
+  //
+  if (gpu::CAMath::Abs(this->getSnp()) >= o2::constants::math::Almost1 ||
+      gpu::CAMath::Abs(this->getTgl()) <= o2::constants::math::Almost0) {
+    return false;
+  }
+  value_T csa, sna, csP, snP, csp = gpu::CAMath::Sqrt((1. - this->getSnp()) * (1. + this->getSnp()));
+  math_utils::detail::sincos(value_T(this->getAlpha()), sna, csa);
+  math_utils::detail::sincos(value_T(t.getPhi()), snP, csP);
+  /*
+    Jacobian is
+    /-sna  -csP/tgL 0     0 0 \
+    | csa  -snP/tgL 0     0 0 |
+    | 0     0       1/csp 0 0 |
+    | 0     0       0     1 0 |
+    \ 0     0       0     0 1 /
+  */
+  auto tgLI = 1 / this->getTgl();
+  const value_T d1 = -sna;
+  const value_T d2 = -csP * tgLI;
+  const value_T e1 = csa;
+  const value_T e2 = -snP * tgLI;
+  const value_T f1 = 1 / csp;
+  SMatrix55Sym C;
+  C(0, 0) = d1 * d1 * getSigmaY2() + 2 * d1 * d2 * getSigmaZY() + d2 * d2 * getSigmaZ2();
+  C(0, 1) = d1 * e1 * getSigmaY2() + (d1 * e2 + d2 * e1) * getSigmaZY() + d2 * e2 * getSigmaZ2();
+  C(1, 1) = e1 * e1 * getSigmaY2() + 2 * e1 * e2 * getSigmaZY() + e2 * e2 * getSigmaZ2();
+
+  C(0, 2) = f1 * (d1 * getSigmaSnpY() + d2 * getSigmaSnpZ());
+  C(1, 2) = f1 * (e1 * getSigmaSnpY() + e2 * getSigmaSnpZ());
+  C(2, 2) = f1 * f1 * getSigmaSnp2();
+
+  C(0, 3) = d1 * getSigmaTglY() + d2 * getSigmaTglZ();
+  C(1, 3) = e1 * getSigmaTglY() + e2 * getSigmaTglZ();
+  C(2, 3) = f1 * getSigmaTglSnp();
+  C(3, 3) = getSigmaTgl2();
+
+  C(0, 4) = d1 * getSigma1PtY() + d2 * getSigma1PtZ();
+  C(1, 4) = e1 * getSigma1PtY() + e2 * getSigma1PtZ();
+  C(2, 4) = f1 * getSigma1PtSnp();
+  C(3, 4) = getSigma1PtTgl();
+  C(4, 4) = getSigma1Pt2();
+  t.setCovariances(C);
+  return true;
+}
+#endif
+
 namespace o2::track
 {
 #if !defined(GPUCA_GPUCODE) || defined(GPUCA_GPUCODE_DEVICE) // FIXME: DR: WORKAROUND to avoid CUDA bug creating host symbols for device code.

Detectors/Calibration/testMacros/getRunParameters.cxx

@@ -53,6 +53,18 @@ void writeDurationToFile(long duration)
   fclose(fptr);
 }
 
+void writeSORToFile(long sor)
+{
+
+  FILE* fptr = fopen("SOR.txt", "w");
+  if (fptr == nullptr) {
+    LOGP(fatal, "ERROR: Could not open file to write SOR!");
+    return;
+  }
+  fprintf(fptr, "%ld", sor);
+  fclose(fptr);
+}
+
 void writeBFieldToFile(float b)
 {
 
@@ -165,6 +177,7 @@ int main(int argc, char* argv[])
       ir = -1.f;
       writeIRtoFile(ir);
       writeDurationToFile(run_O2duration);
+      writeSORToFile(tsSOR);
       return 0;
     }
   }
@@ -204,6 +217,7 @@ int main(int argc, char* argv[])
   }
   writeIRtoFile(ir);
   writeDurationToFile(duration);
+  writeSORToFile(tsSOR);
 
   return 0;
 }

Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu

@@ -146,7 +146,7 @@ GPUdii() o2::track::TrackParCov buildTrackSeed(const Cluster& cluster1,
     const float dx = x3 - x1;
     const float dy = y3 - y1;
     snp = sign * dy / o2::gpu::CAMath::Hypot(dx, dy);
-    q2pt = sign / track::kMostProbablePt;
+    q2pt = 1.f / track::kMostProbablePt;
     q2pt2 = 1.f;
   } else {
     const float crv = math_utils::computeCurvature(x3, y3, x2, y2, x1, y1);

Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx

@@ -1281,7 +1281,7 @@ track::TrackParCov TrackerTraits<nLayers>::buildTrackSeed(const Cluster& cluster
     const float dx = x3 - x1;
     const float dy = y3 - y1;
     snp = sign * dy / o2::gpu::CAMath::Hypot(dx, dy);
-    q2pt = sign / track::kMostProbablePt;
+    q2pt = 1.f / track::kMostProbablePt;
     q2pt2 = 1.f;
   } else {
     const float crv = math_utils::computeCurvature(x3, y3, x2, y2, x1, y1);