Merge branch 'dev' into onnx_gpu_timer

Author: ChSonnabend

DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackFwd.h

@@ -161,6 +161,7 @@ class TrackParCovFwd : public TrackParFwd
   void propagateToZquadratic(double zEnd, double zField);
   void propagateToZhelix(double zEnd, double zField);
   void propagateToZ(double zEnd, double zField); // Parameters: helix; errors: quadratic
+  void propagateToDCAhelix(double zField, const std::array<double, 3>& p, std::array<double, 3>& dca);
 
   // Add Multiple Coulomb Scattering effects
   void addMCSEffect(double x2X0);

DataFormats/Reconstruction/include/ReconstructionDataFormats/Vertex.h

@@ -18,10 +18,13 @@
 
 #include "CommonDataFormat/TimeStamp.h"
 #ifndef GPUCA_GPUCODE_DEVICE
-#include <iosfwd>
-#include <string>
 #include <type_traits>
 #include <array>
+#ifndef GPUCA_NO_FMT
+#include <sstream>
+#include <string>
+#include <fmt/format.h>
+#endif
 #endif
 
 namespace o2
@@ -135,6 +138,11 @@ class Vertex : public VertexBase
   {
   }
 
+#if !defined(GPUCA_NO_FMT) && !defined(GPUCA_GPUCODE_DEVICE)
+  void print() const;
+  std::string asString() const;
+#endif
+
   GPUd() ushort getNContributors() const { return mNContributors; }
   GPUd() void setNContributors(ushort v) { mNContributors = v; }
   GPUd() void addContributor() { mNContributors++; }
@@ -162,6 +170,49 @@ class Vertex : public VertexBase
 
 #if !defined(GPUCA_GPUCODE_DEVICE) && !defined(GPUCA_NO_FMT)
 std::ostream& operator<<(std::ostream& os, const o2::dataformats::VertexBase& v);
+
+namespace detail
+{
+template <typename T>
+concept Streamable = requires(std::ostream& os, const T& a) {
+  { os << a } -> std::same_as<std::ostream&>;
+};
+
+template <typename T>
+concept HasFormattableTimeStamp = requires(const T& t) {
+  { fmt::format("{}", t.getTimeStamp()) } -> std::convertible_to<std::string>;
+};
+} // namespace detail
+
+template <typename Stamp>
+inline std::string Vertex<Stamp>::asString() const
+{
+  const std::string stamp = [&]() -> std::string {
+    if constexpr (detail::Streamable<Stamp>) {
+      std::ostringstream oss;
+      oss << mTimeStamp;
+      return oss.str();
+    } else if constexpr (detail::HasFormattableTimeStamp<Stamp>) {
+      return fmt::format("{}", mTimeStamp.getTimeStamp());
+    } else {
+      return "X";
+    }
+  }();
+  return fmt::format("{} NContrib:{} Chi2:{:.2f} Flags:{:b} Stamp:{}", VertexBase::asString(), mNContributors, mChi2, mBits, stamp);
+}
+
+template <typename Stamp>
+inline std::ostream& operator<<(std::ostream& os, const o2::dataformats::Vertex<Stamp>& v)
+{
+  os << v.asString();
+  return os;
+}
+
+template <typename Stamp>
+inline void Vertex<Stamp>::print() const
+{
+  std::cout << *this << '\n';
+}
 #endif
 
 } // namespace dataformats

DataFormats/Reconstruction/src/TrackFwd.cxx

@@ -11,6 +11,7 @@
 
 #include "ReconstructionDataFormats/TrackFwd.h"
 #include "Math/MatrixFunctions.h"
+#include <GPUCommonLogger.h>
 
 namespace o2
 {
@@ -503,5 +504,73 @@ bool TrackParCovFwd::getCovXYZPxPyPzGlo(std::array<float, 21>& cv) const
   return true;
 }
 
+//________________________________________________________________
+
+void TrackParCovFwd::propagateToDCAhelix(double zField, const std::array<double, 3>& p, std::array<double, 3>& dca)
+{
+  // Computing DCA of fwd track w.r.t vertex in helix track model, using Newton-Raphson minimization
+
+  auto x0 = mParameters(0);
+  auto y0 = mParameters(1);
+  auto z0 = mZ;
+  auto phi0 = mParameters(2);
+  auto tanl = mParameters(3);
+  auto qOverPt = mParameters(4);
+  auto k = TMath::Abs(o2::constants::math::B2C * zField);
+  auto qpt = 1.0 / qOverPt;
+  auto qR = qpt / std::fabs(k);
+  auto invtanl = 1.0 / tanl;
+  auto Hz = std::copysign(1, zField);
+
+  auto xPV = p[0];
+  auto yPV = p[1];
+  auto zPV = p[2];
+
+  auto qRtanl = qR * tanl;
+  auto invqRtanl = 1.0 / qRtanl;
+  auto [sinp, cosp] = o2::math_utils::sincosd(phi0);
+
+  auto z = zPV;
+  double tol = 1e-4;
+  int max_iter = 10;
+  int iter = 0;
+
+  while (iter++ < max_iter) {
+    double theta = (z0 - z) * invqRtanl;
+    double phi_theta = phi0 + Hz * theta;
+    double sin_phi_theta = sin(phi_theta);
+    double cos_phi_theta = cos(phi_theta);
+
+    double DX = x0 - Hz * qR * (sin_phi_theta - sinp) - xPV;
+    double DY = y0 + Hz * qR * (cos_phi_theta - cosp) - yPV;
+    double DZ = z - zPV;
+
+    double dD2_dZ =
+      2 * DX * cos_phi_theta * invtanl +
+      2 * DY * sin_phi_theta * invtanl +
+      2 * DZ;
+
+    double d2D2_dZ2 =
+      2 * invtanl * invtanl +
+      2 * invtanl * (DX * Hz * sin_phi_theta - DY * Hz * cos_phi_theta) * invqRtanl +
+      2;
+
+    double z_new = z - dD2_dZ / d2D2_dZ2;
+
+    if (std::abs(z_new - z) < tol) {
+      z = z_new;
+      this->propagateToZhelix(z, zField);
+      dca[0] = this->getX() - xPV;
+      dca[1] = this->getY() - yPV;
+      dca[2] = this->getZ() - zPV;
+      LOG(debug) << "Converged after " << iter << " iterations for vertex X=" << p[0] << ", Y=" << p[1] << ", Z = " << p[2];
+      return;
+    }
+    z = z_new;
+  }
+  LOG(debug) << "Failed to converge after " << iter << " iterations for vertex X=" << p[0] << ", Y=" << p[1] << ", Z = " << p[2];
+  return;
+}
+
 } // namespace track
 } // namespace o2

DataFormats/Reconstruction/src/Vertex.cxx

@@ -10,9 +10,9 @@
 // or submit itself to any jurisdiction.
 
 #include "ReconstructionDataFormats/Vertex.h"
-#include <iostream>
 #ifndef GPUCA_NO_FMT
-#include <fmt/printf.h>
+#include <iostream>
+#include <fmt/format.h>
 #endif
 
 namespace o2

Detectors/TPC/calibration/include/TPCCalibration/CorrMapParam.h

@@ -29,6 +29,7 @@ struct CorrMapParam : public o2::conf::ConfigurableParamHelper<CorrMapParam> {
   float lumiMean = 0.;       // override TPC corr.map mean lumi (if > 0), disable corrections if < 0
   float lumiMeanRef = 0.;    // override TPC corr.mapRef mean lumi (if > 0)"
   float lumiInstFactor = 1.; // scaling to apply to instantaneous lumi from CTP (but not to IDC scaler)
+  float CTP2IDCFallBackThreshold = 30.; // if needed, interpret map->getLumi() as map->getIDC(), provided map->getLumi() is below this threshold
   int ctpLumiSource = 0;     // CTP lumi source: 0 = LumiInfo.getLumi(), 1 = LumiInfo.getLumiAlt()
 
   O2ParamDef(CorrMapParam, "TPCCorrMap");

Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h

@@ -63,6 +63,7 @@ class CorrectionMapsLoader : public o2::gpu::CorrectionMapsHelper
   void init(o2::framework::InitContext& ic);
   void copySettings(const CorrectionMapsLoader& src);
   void updateInverse(); /// recalculate inverse correction
+  float getMapMeanRate(const o2::gpu::TPCFastTransform* mp, bool lumiOverridden) const;
 
   static void requestCCDBInputs(std::vector<o2::framework::InputSpec>& inputs, std::vector<o2::framework::ConfigParamSpec>& options, const CorrectionMapsLoaderGloOpts& gloOpts);
   static void addGlobalOptions(std::vector<o2::framework::ConfigParamSpec>& options);

Detectors/TPC/calibration/src/CorrectdEdxDistortions.cxx

@@ -70,8 +70,8 @@ void o2::tpc::CorrectdEdxDistortions::setLumi(float lumi)
     LOGP(warn, "Nullptr detected in accessing the correction maps");
     return;
   }
-  const float lumiAvg = mCorrAvg->getLumi();
-  const float lumiDer = mCorrDer->getLumi();
+  const float lumiAvg = mCorrAvg->getIDC();
+  const float lumiDer = mCorrDer->getIDC();
   mScaleDer = (lumi - lumiAvg) / lumiDer;
   LOGP(info, "Setting mScaleDer: {} for inst lumi: {}  avg lumi: {}  deriv. lumi: {}", mScaleDer, lumi, lumiAvg, lumiDer);
 }

Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx

@@ -176,20 +176,52 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher,
   if (matcher == ConcreteDataMatcher("TPC", "CorrMap", 0)) {
     setCorrMap((o2::gpu::TPCFastTransform*)obj);
     mCorrMap->rectifyAfterReadingFromFile();
-    if (getMeanLumiOverride() == 0 && mCorrMap->getLumi() > 0.) {
-      setMeanLumi(mCorrMap->getLumi(), false);
+    mCorrMap->setCTP2IDCFallBackThreshold(o2::tpc::CorrMapParam::Instance().CTP2IDCFallBackThreshold);
+    if (getMeanLumiOverride() != 0) {
+      if (getLumiScaleType() == 1) {
+        mCorrMap->setLumi(getMeanLumiOverride());
+        LOGP(info, "CorrMap mean lumi rate is overridden to {}", mCorrMap->getLumi());
+      } else if (getLumiScaleType() == 2) {
+        mCorrMap->setIDC(getMeanLumiOverride());
+        LOGP(info, "CorrMap mean IDC rate is overridden to {}", mCorrMap->getIDC());
+      }
+    }
+    float mapMeanRate = 0;
+    if (getLumiScaleType() == 1) {
+      mapMeanRate = mCorrMap->getLumi();
+    } else if (getLumiScaleType() == 2) {
+      mapMeanRate = mCorrMap->getIDC();
+    }
+    if (getMeanLumiOverride() == 0 && mapMeanRate > 0.) {
+      setMeanLumi(mapMeanRate, false);
     }
-    LOGP(debug, "MeanLumiOverride={} MeanLumiMap={} -> meanLumi = {}", getMeanLumiOverride(), mCorrMap->getLumi(), getMeanLumi());
+    LOGP(debug, "MeanLumiOverride={} MeanLumiMap={} -> meanLumi = {}", getMeanLumiOverride(), mapMeanRate, getMeanLumi());
     setUpdatedMap();
     return true;
   }
   if (matcher == ConcreteDataMatcher("TPC", "CorrMapRef", 0)) {
     setCorrMapRef((o2::gpu::TPCFastTransform*)obj);
     mCorrMapRef->rectifyAfterReadingFromFile();
+    mCorrMapRef->setCTP2IDCFallBackThreshold(o2::tpc::CorrMapParam::Instance().CTP2IDCFallBackThreshold);
+    if (getMeanLumiRefOverride() != 0) {
+      if (getLumiScaleType() == 1) {
+        mCorrMapRef->setLumi(getMeanLumiRefOverride());
+        LOGP(info, "CorrMapRef mean lumi rate is overridden to {}", mCorrMapRef->getLumi());
+      } else if (getLumiScaleType() == 2) {
+        mCorrMapRef->setIDC(getMeanLumiRefOverride());
+        LOGP(info, "CorrMapRef mean IDC rate is overridden to {}", mCorrMapRef->getIDC());
+      }
+    }
+    float mapRefMeanRate = 0;
+    if (getLumiScaleType() == 1) {
+      mapRefMeanRate = mCorrMapRef->getLumi();
+    } else if (getLumiScaleType() == 2) {
+      mapRefMeanRate = mCorrMapRef->getIDC();
+    }
     if (getMeanLumiRefOverride() == 0) {
-      setMeanLumiRef(mCorrMapRef->getLumi());
+      setMeanLumiRef(mapRefMeanRate);
     }
-    LOGP(debug, "MeanLumiRefOverride={} MeanLumiMap={} -> meanLumi = {}", getMeanLumiRefOverride(), mCorrMapRef->getLumi(), getMeanLumiRef());
+    LOGP(debug, "MeanLumiRefOverride={} MeanLumiMap={} -> meanLumi = {}", getMeanLumiRefOverride(), mapRefMeanRate, getMeanLumiRef());
     setUpdatedMapRef();
     return true;
   }
@@ -217,7 +249,7 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher,
     int scaleType = getLumiScaleType();
     const std::array<std::string, 3> lumiS{"OFF", "CTP", "TPC scaler"};
     if (scaleType >= lumiS.size()) {
-      LOGP(fatal, "Wrong lumi-scale-type provided!");
+      LOGP(fatal, "Wrong corrmap-lumi-mode provided!");
     }
 
     LOGP(info, "TPC correction map params updated: SP corrections: {} (corr.map scaling type={}, override values: lumiMean={} lumiRefMean={} lumiScaleMode={}), CTP Lumi: source={} lumiInstOverride={} , LumiInst scale={} ",

Detectors/TRD/workflow/include/TRDWorkflow/GainCalibSpec.h

@@ -45,7 +45,7 @@ class GainCalibDevice : public o2::framework::Task
   {
     o2::base::GRPGeomHelper::instance().setRequest(mCCDBRequest);
     auto slotL = ic.options().get<uint32_t>("sec-per-slot");
-    auto delay = ic.options().get<uint32_t>("max-delay");
+    auto delay = ic.options().get<float>("max-delay");
     mCalibrator = std::make_unique<o2::trd::CalibratorGain>();
     mCalibrator->setSlotLengthInSeconds(slotL);
     mCalibrator->setMaxSlotsDelay(delay);
@@ -155,7 +155,7 @@ DataProcessorSpec getTRDGainCalibSpec()
     AlgorithmSpec{adaptFromTask<o2::calibration::GainCalibDevice>(ccdbRequest)},
     Options{
       {"sec-per-slot", VariantType::UInt32, 900u, {"number of seconds per calibration time slot"}},
-      {"max-delay", VariantType::UInt32, 2u, {"number of slots in past to consider"}},
+      {"max-delay", VariantType::Float, 0.05f, {"number of slots in past to consider"}},
       {"enable-root-output", VariantType::Bool, false, {"output tprofiles and fits to root file"}},
     }};
 }

Detectors/TRD/workflow/include/TRDWorkflow/T0FitSpec.h

@@ -46,7 +46,7 @@ class T0FitDevice : public o2::framework::Task
   {
     o2::base::GRPGeomHelper::instance().setRequest(mCCDBRequest);
     auto slotL = ic.options().get<uint32_t>("sec-per-slot");
-    auto delay = ic.options().get<uint32_t>("max-delay");
+    auto delay = ic.options().get<float>("max-delay");
 
     mFitInstance = std::make_unique<o2::trd::T0Fit>();
     mFitInstance->setSlotLengthInSeconds(slotL);
@@ -159,7 +159,7 @@ DataProcessorSpec getTRDT0FitSpec()
     AlgorithmSpec{adaptFromTask<device>(ccdbRequest)},
     Options{
       {"sec-per-slot", VariantType::UInt32, 900u, {"number of seconds per calibration time slot"}},
-      {"max-delay", VariantType::UInt32, 2u, {"number of slots in past to consider"}},
+      {"max-delay", VariantType::Float, 0.05f, {"number of slots in past to consider"}},
       {"enable-root-output", VariantType::Bool, false, {"output t0 values to root file"}},
     }};
 }