Merge branch 'dev' into its/capcatiy

Author: f3sch

Common/DCAFitter/GPU/cuda/CMakeLists.txt

@@ -20,6 +20,9 @@ o2_add_library(DCAFitterCUDA
                                      O2::DetectorsBase
                PRIVATE_LINK_LIBRARIES O2::GPUTrackingCUDAExternalProvider)
 set_property(TARGET ${targetName} PROPERTY CUDA_SEPARABLE_COMPILATION ON)
+# Device LTO, so that the device link can inline across the
+# O2::GPUTrackingCUDAExternalProvider objects, which are compiled to LTO IR.
+set_property(TARGET ${targetName} PROPERTY INTERPROCEDURAL_OPTIMIZATION ON)
 # add_compile_options(-lineinfo)
 
 #o2_add_test(DCAFitterNCUDA

Common/Field/include/Field/MagneticField.h

@@ -262,8 +262,8 @@ class MagneticField : public FairField
 
   Int_t mDefaultIntegration;             ///< Default integration method as indicated in Geant
   Int_t mPrecisionInteg;                 ///< Alternative integration method, e.g. for higher precision
-  Double_t mMultipicativeFactorSolenoid; ///< Multiplicative factor for solenoid
-  Double_t mMultipicativeFactorDipole;   ///< Multiplicative factor for dipole
+  Double_t mMultipicativeFactorSolenoid; ///< Multiplicative factor for solenoid, polarity convention applied
+  Double_t mMultipicativeFactorDipole;   ///< Multiplicative factor for dipole, polarity convention applied
   Double_t mMaxField;                    ///< Max Field as indicated in Geant
   Bool_t mDipoleOnOffFlag;               ///< Dipole ON/OFF flag
 

Common/Field/src/MagneticField.cxx

@@ -108,8 +108,8 @@ MagneticField::MagneticField(const char* name, const char* title, Double_t facto
     mBeamEnergy(be),
     mDefaultIntegration(integ),
     mPrecisionInteg(1),
-    mMultipicativeFactorSolenoid(factorSol),
-    mMultipicativeFactorDipole(factorDip),
+    mMultipicativeFactorSolenoid(1.),
+    mMultipicativeFactorDipole(1.),
     mMaxField(fmax),
     mDipoleOnOffFlag(factorDip == 0.),
     mQuadrupoleGradient(0),
@@ -122,6 +122,8 @@ MagneticField::MagneticField(const char* name, const char* title, Double_t facto
   /*
    * Constructor for human readable params
    */
+  setFactorSolenoid(factorSol);
+  setFactorDipole(factorDip);
   setDataFileName(path.c_str());
   if (!gOriginBias) {
     checkOriginBias();
@@ -139,8 +141,8 @@ MagneticField::MagneticField(const MagFieldParam& param)
     mBeamEnergy(param.GetBeamEnergy()),
     mDefaultIntegration(param.GetDefInt()),
     mPrecisionInteg(1),
-    mMultipicativeFactorSolenoid(param.GetFactorSol()), // temporary
-    mMultipicativeFactorDipole(param.GetFactorDip()),   // temporary
+    mMultipicativeFactorSolenoid(1.),
+    mMultipicativeFactorDipole(1.),
     mMaxField(param.GetMaxField()),
     mDipoleOnOffFlag(param.GetFactorDip() == 0.),
     mQuadrupoleGradient(0),
@@ -153,6 +155,8 @@ MagneticField::MagneticField(const MagFieldParam& param)
   /*
    * Constructor for FairParam derived params
    */
+  setFactorSolenoid(param.GetFactorSol());
+  setFactorDipole(param.GetFactorDip());
   setDataFileName(param.GetMapPath());
   if (!gOriginBias) {
     checkOriginBias();
@@ -235,8 +239,11 @@ void MagneticField::CreateField()
 
   loadParameterization();
   initializeMachineField(mBeamType, mBeamEnergy);
-  setFactorSolenoid(mMultipicativeFactorSolenoid);
-  setFactorDipole(mMultipicativeFactorDipole);
+  // The scaling factors are left alone: they already carry the polarity convention, and
+  // re-applying it would invert a field re-initialized after being read back from a file.
+  if (mFastField) {
+    mFastField->setFactorSol(getFactorSolenoid());
+  }
   double xyz[3] = {0., 0., 0.};
   mSolenoid = getBz(xyz);
   Print("a");

Common/Field/test/testMagneticField.cxx

@@ -18,6 +18,7 @@
 #include "Field/MagFieldFast.h"
 #include <memory>
 #include <fairlogger/Logger.h> // for FairLogger
+#include <TFile.h>
 #include <TStopwatch.h>
 #include <TRandom.h>
 
@@ -98,3 +99,48 @@ BOOST_AUTO_TEST_CASE(MagneticField_test)
     BOOST_CHECK(TMath::Abs(rms[i] / nomBz) < 1.e-3);
   }
 }
+
+BOOST_AUTO_TEST_CASE(MagneticField_reinitialization_test)
+{
+  // The measured map is transient, so a MagneticField read back from a file has to be
+  // re-created before it can be used. That must reproduce the field vectors, not merely
+  // their magnitude: a sign flip leaves |B| untouched.
+  const double points[][3] = {
+    {0., 0., 0.},      // solenoid, on axis
+    {100., 50., 100.}, // solenoid, off axis
+    {10., 10., -900.}, // muon dipole
+    {0., 0., 1000.},   // compensator 1A, side A
+    {0., 0., -2049.},  // compensator 2C, side C
+    {0., 0., 2049.}    // compensator 2A, side A
+  };
+  const int npoints = sizeof(points) / sizeof(points[0]);
+  const double tolerance = 1.e-9; // kGauss
+
+  std::unique_ptr<MagneticField> fld = std::make_unique<MagneticField>("Maps", "Maps", 1., 1., MagFieldParam::k5kG);
+  const double facSol = fld->getFactorSolenoid(), facDip = fld->getFactorDipole();
+  double bref[npoints][3] = {};
+  for (int ip = 0; ip < npoints; ip++) {
+    fld->Field(points[ip], bref[ip]);
+    // a point where the field vanishes would make the comparisons below vacuous
+    BOOST_CHECK(TMath::Abs(bref[ip][0]) + TMath::Abs(bref[ip][1]) + TMath::Abs(bref[ip][2]) > tolerance);
+  }
+
+  const char* fname = "testMagneticFieldReinitialization.root";
+  {
+    TFile fout(fname, "recreate");
+    fout.WriteObject(fld.get(), "field");
+  }
+  TFile fin(fname);
+  auto* fldRead = fin.Get<MagneticField>("field");
+  BOOST_REQUIRE(fldRead != nullptr);
+  fldRead->CreateField();
+  BOOST_CHECK_EQUAL(fldRead->getFactorSolenoid(), facSol);
+  BOOST_CHECK_EQUAL(fldRead->getFactorDipole(), facDip);
+  for (int ip = 0; ip < npoints; ip++) {
+    double b[3] = {};
+    fldRead->Field(points[ip], b);
+    for (int i = 0; i < 3; i++) {
+      BOOST_CHECK_SMALL(b[i] - bref[ip][i], tolerance);
+    }
+  }
+}

Common/MathUtils/include/MathUtils/Cartesian.h

@@ -38,6 +38,7 @@
 
 #include "GPUROOTCartesianFwd.h"
 #include "GPUROOTSMatrixFwd.h"
+#include "MathUtils/detail/trigonometric.h"
 
 namespace o2
 {
@@ -68,7 +69,7 @@ class Rotation2D
 
   Rotation2D() = default;
   Rotation2D(value_t cs, value_t sn) : mCos(cs), mSin(sn) {}
-  Rotation2D(value_t phiZ) : mCos(cos(phiZ)), mSin(sin(phiZ)) {}
+  Rotation2D(value_t phiZ) { detail::sincos<value_t>(phiZ, mSin, mCos); }
   ~Rotation2D() = default;
   Rotation2D(const Rotation2D& src) = default;
   Rotation2D(Rotation2D&& src) = default;
@@ -77,8 +78,7 @@ class Rotation2D
 
   void set(value_t phiZ)
   {
-    mCos = cos(phiZ);
-    mSin = sin(phiZ);
+    detail::sincos<value_t>(phiZ, mSin, mCos);
   }
 
   void set(value_t cs, value_t sn)

Common/SimConfig/include/SimConfig/SimParams.h

@@ -23,6 +23,7 @@ namespace conf
 // (mostly used in O2MCApplication stepping)
 struct SimCutParams : public o2::conf::ConfigurableParamHelper<SimCutParams> {
   bool stepFiltering = true; // if we activate the step filtering in O2BaseMCApplication
+  std::string stepFilteringMacro = "";                                                                        // ROOT macro providing keepStep(); empty = built-in z/R cut
   bool stepTrackRefHook = false;                                                                              // if we create track references during generic stepping
   std::string stepTrackRefHookFile = "${O2_ROOT}/share/Detectors/gconfig/StandardSteppingTrackRefHook.macro"; // the standard code holding the TrackRef callback
 

DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h

@@ -548,7 +548,7 @@ GPUdi() void TrackParametrization<value_T>::getLineParams(o2::math_utils::Interv
 template <typename value_T>
 GPUdi() auto TrackParametrization<value_T>::getCurvature(value_t b) const -> value_t
 {
-  return mAbsCharge ? mP[kQ2Pt] * b * o2::constants::math::B2C : 0.;
+  return mAbsCharge ? mP[kQ2Pt] * b * o2::constants::math::B2C : value_T(0);
 }
 
 //____________________________________________________________

DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h

@@ -114,6 +114,7 @@ class TrackParametrizationWithError : public TrackParametrization<value_T>
   GPUd() void buildCombinedCovMatrix(const TrackParametrizationWithError& rhs, MatrixDSym5& cov) const;
   GPUd() value_t getPredictedChi2(const TrackParametrizationWithError& rhs, MatrixDSym5& covToSet) const;
   GPUd() value_t getPredictedChi2(const TrackParametrizationWithError& rhs) const;
+  GPUd() value_t getPredictedChi2Fast(const TrackParametrizationWithError& rhs) const;
   GPUd() value_t getPredictedChi2Quiet(const TrackParametrizationWithError& rhs) const;
   GPUd() bool update(const TrackParametrizationWithError& rhs, const MatrixDSym5& covInv);
   GPUd() bool update(const TrackParametrizationWithError& rhs);

DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h

@@ -66,26 +66,26 @@ GPUd() void g3helx3(value_T qfield, value_T step, std::array<value_T, 7>& vect)
   static_assert(std::is_floating_point_v<value_T>);
 #endif
 
-  const int ix = 0, iy = 1, iz = 2, ipx = 3, ipy = 4, ipz = 5, ipp = 6;
-  constexpr value_T kOvSqSix = 0.408248f; // std::sqrt(1./6.);
+  constexpr int ix = 0, iy = 1, iz = 2, ipx = 3, ipy = 4, ipz = 5, ipp = 6;
+  constexpr value_T kOvSqSix = value_T(0.408248); // std::sqrt(1./6.);
 
   value_T cosx = vect[ipx], cosy = vect[ipy], cosz = vect[ipz];
 
   value_T rho = qfield * constants::math::B2C / vect[ipp];
   value_T tet = rho * step;
 
   value_T tsint, sintt, sint, cos1t;
-  if (gpu::CAMath::Abs(tet) > 0.03f) {
+  if (gpu::CAMath::Abs(tet) > value_T(0.03)) {
     sint = gpu::CAMath::Sin(tet);
     sintt = sint / tet;
     tsint = (tet - sint) / tet;
-    value_T t = gpu::CAMath::Sin(0.5f * tet);
-    cos1t = 2 * t * t / tet;
+    value_T t = gpu::CAMath::Sin(value_T(0.5) * tet);
+    cos1t = value_T(2) * t * t / tet;
   } else {
-    tsint = tet * tet / 6.f;
-    sintt = (1.f - tet * kOvSqSix) * (1.f + tet * kOvSqSix); // 1.- tsint;
+    tsint = tet * tet / value_T(6);
+    sintt = (value_T(1) - tet * kOvSqSix) * (value_T(1) + tet * kOvSqSix); // 1.- tsint;
     sint = tet * sintt;
-    cos1t = 0.5f * tet;
+    cos1t = value_T(0.5) * tet;
   }
 
   value_T f1 = step * sintt;
@@ -124,25 +124,25 @@ GPUd() value_T BetheBlochSolid(value_T bg, value_T rho, value_T kp1, value_T kp2
   static_assert(std::is_floating_point_v<value_T>);
 #endif
 
-  constexpr value_T mK = 0.307075e-3; // [GeV*cm^2/g]
-  constexpr value_T me = 0.511e-3;    // [GeV/c^2]
-  kp1 *= 2.303f;
-  kp2 *= 2.303f;
-  value_T bg2 = bg * bg, beta2 = bg2 / (1 + bg2);
-  value_T maxT = 2.f * me * bg2; // neglecting the electron mass
+  constexpr value_T mK = value_T(0.307075e-3); // [GeV*cm^2/g]
+  constexpr value_T me = value_T(0.511e-3);    // [GeV/c^2]
+  kp1 *= value_T(2.303);
+  kp2 *= value_T(2.303);
+  value_T bg2 = bg * bg, beta2 = bg2 / (value_T(1) + bg2);
+  value_T maxT = value_T(2) * me * bg2; // neglecting the electron mass
 
   //*** Density effect
-  value_T d2 = 0.;
+  value_T d2 = value_T(0);
   const value_T x = gpu::CAMath::Log(bg);
-  const value_T lhwI = gpu::CAMath::Log(28.816f * 1e-9f * gpu::CAMath::Sqrt(rho * meanZA) / meanI);
+  const value_T lhwI = gpu::CAMath::Log(value_T(28.816e-9) * gpu::CAMath::Sqrt(rho * meanZA) / meanI);
   if (x > kp2) {
-    d2 = lhwI + x - 0.5f;
+    d2 = lhwI + x - value_T(0.5);
   } else if (x > kp1) {
     double r = (kp2 - x) / (kp2 - kp1);
-    d2 = lhwI + x - 0.5f + (0.5f - lhwI - kp1) * r * r * r;
+    d2 = lhwI + x - value_T(0.5) + (value_T(0.5) - lhwI - kp1) * r * r * r;
   }
-  auto dedx = mK * meanZA / beta2 * (0.5f * gpu::CAMath::Log(2 * me * bg2 * maxT / (meanI * meanI)) - beta2 - d2);
-  return dedx > 0. ? dedx : 0.;
+  auto dedx = mK * meanZA / beta2 * (value_T(0.5) * gpu::CAMath::Log(value_T(2) * me * bg2 * maxT / (meanI * meanI)) - beta2 - d2);
+  return dedx > value_T(0) ? dedx : value_T(0);
 }
 
 //____________________________________________________
@@ -166,26 +166,26 @@ GPUd() value_T BetheBlochSolidOpt(value_T bg)
   //  constexpr value_T meanI = 173e-9;
   //  constexpr value_T me = 0.511e-3;    // [GeV/c^2]
 
-  constexpr value_T mK = 0.307075e-3; // [GeV*cm^2/g]
-  constexpr value_T kp1 = 0.20 * 2.303;
-  constexpr value_T kp2 = 3.00 * 2.303;
-  constexpr value_T meanZA = 0.49848;
-  constexpr value_T lhwI = -1.7175226;         // gpu::CAMath::Log(28.816 * 1e-9 * gpu::CAMath::Sqrt(rho * meanZA) / meanI);
-  constexpr value_T log2muTomeanI = 8.6839805; // gpu::CAMath::Log( 2. * me / meanI);
+  constexpr value_T mK = value_T(0.307075e-3); // [GeV*cm^2/g]
+  constexpr value_T kp1 = value_T(0.20 * 2.303);
+  constexpr value_T kp2 = value_T(3.00 * 2.303);
+  constexpr value_T meanZA = value_T(0.49848);
+  constexpr value_T lhwI = value_T(-1.7175226);         // gpu::CAMath::Log(28.816 * 1e-9 * gpu::CAMath::Sqrt(rho * meanZA) / meanI);
+  constexpr value_T log2muTomeanI = value_T(8.6839805); // gpu::CAMath::Log( 2. * me / meanI);
 
-  value_T bg2 = bg * bg, beta2 = bg2 / (1. + bg2);
+  value_T bg2 = bg * bg, beta2 = bg2 / (value_T(1) + bg2);
 
   //*** Density effect
-  value_T d2 = 0.;
+  value_T d2 = value_T(0);
   const value_T x = gpu::CAMath::Log(bg);
   if (x > kp2) {
-    d2 = lhwI - 0.5f + x;
+    d2 = lhwI - value_T(0.5) + x;
   } else if (x > kp1) {
     value_T r = (kp2 - x) / (kp2 - kp1);
-    d2 = lhwI - 0.5 + x + (0.5 - lhwI - kp1) * r * r * r;
+    d2 = lhwI - value_T(0.5) + x + (value_T(0.5) - lhwI - kp1) * r * r * r;
   }
   auto dedx = mK * meanZA / beta2 * (log2muTomeanI + x + x - beta2 - d2);
-  return dedx > 0. ? dedx : 0.;
+  return dedx > value_T(0) ? dedx : value_T(0);
 }
 
 //____________________________________________________
@@ -203,12 +203,12 @@ GPUdi() value_T BetheBlochSolidDerivative(value_T dedx, value_T bg)
   // dedx - precalculate dedx for bg
   // bg  - beta*gamma
   //
-  constexpr value_T mK = 0.307075e-3; // [GeV*cm^2/g]
-  constexpr value_T meanZA = 0.49848;
+  constexpr value_T mK = value_T(0.307075e-3); // [GeV*cm^2/g]
+  constexpr value_T meanZA = value_T(0.49848);
   auto bg2 = bg * bg;
-  auto t1 = 1 + bg2;
+  auto t1 = value_T(1) + bg2;
   //  auto derH = (mK * meanZA * (t1+bg2) - dedx*bg2)/(bg*t1);
-  auto derH = (mK * meanZA * (t1 + 1. / bg2) - dedx) / (bg * t1);
+  auto derH = (mK * meanZA * (t1 + value_T(1) / bg2) - dedx) / (bg * t1);
   return derH + derH;
 }
 

DataFormats/Reconstruction/src/TrackParametrization.cxx

@@ -364,7 +364,7 @@ GPUd() bool TrackParametrization<value_T>::propagateParamTo(value_t xk, value_t
   if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
     return false;
   }
-  double dy2dx = (f1 + f2) / (r1 + r2);
+  value_t dy2dx = (f1 + f2) / (r1 + r2);
   bool arcz = gpu::CAMath::Abs(x2r) > 0.05f;
   if (arcz) {
     // for small dx/R the linear apporximation of the arc by the segment is OK,