Merge c389c71 into sapling-pr-archive-ktf

Author: ktf

CODEOWNERS

@@ -28,9 +28,9 @@
 /DataFormats/Detectors/Common                  @shahor02
 /DataFormats/Detectors/CPV                     @peressounko @kharlov
 /DataFormats/Detectors/CTP                     @lietava
-/DataFormats/Detectors/EMCAL                   @mfasDa @jokonig
+/DataFormats/Detectors/EMCAL                   @nstrangm @jokonig
 /DataFormats/Detectors/FIT                     @jotwinow @afurs @andreasmolander @sahilupadhyaya92
-/DataFormats/Detectors/FOCAL                   @maxrauch @mfasDa @iarsene @matthiasrichter
+/DataFormats/Detectors/FOCAL                   @maxrauch @iarsene @matthiasrichter
 /DataFormats/Detectors/GlobalTracking          @shahor02
 /DataFormats/Detectors/GlobalTrackingWorkflow  @shahor02
 /DataFormats/Detectors/HMPID                   @gvolpe79
@@ -58,9 +58,9 @@
 /Detectors/Base                    @sawenzel @shahor02
 /Detectors/Calibration             @chiarazampolli @shahor02
 /Detectors/CPV                     @peressounko @kharlov
-/Detectors/EMCAL                   @mfasDa @jokonig
+/Detectors/EMCAL                   @nstrangm @jokonig
 /Detectors/FIT                     @jotwinow @afurs @andreasmolander @sahilupadhyaya92
-/Detectors/FOCAL                   @maxrauch @mfasDa @iarsene @matthiasrichter
+/Detectors/FOCAL                   @maxrauch @iarsene @matthiasrichter
 /Detectors/Geometry                @sawenzel @shahor02
 /Detectors/GlobalTracking          @shahor02
 /Detectors/GlobalTrackingWorkflow  @shahor02

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/EMCAL/calibration/CMakeLists.txt

@@ -20,6 +20,7 @@ o2_add_library(EMCALCalibration
                         src/PedestalCalibDevice.cxx
                         src/PedestalProcessorDevice.cxx
                         src/PedestalProcessorData.cxx
+                        src/EMCALTempCalibExtractor.cxx
                PUBLIC_LINK_LIBRARIES O2::CCDB O2::EMCALBase
                                      O2::EMCALCalib
                                      O2::EMCALReconstruction
@@ -46,6 +47,7 @@ o2_target_root_dictionary(EMCALCalibration
                                   include/EMCALCalibration/EMCDCSProcessor.h
                                   include/EMCALCalibration/EMCALPedestalHelper.h
                                   include/EMCALCalibration/PedestalProcessorData.h
+                                  include/EMCALCalibration/EMCALTempCalibExtractor.h
                           LINKDEF src/EMCALCalibrationLinkDef.h)
 
 o2_add_executable(emcal-channel-calib-workflow

Detectors/EMCAL/calibration/include/EMCALCalibration/EMCALTempCalibExtractor.h

@@ -0,0 +1,93 @@
+// 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.
+
+/// \class EMCALTempCalibExtractor
+/// \brief  Calculate gain correction factors based on the temperature and the cell-by-cell temperature dependent correction factors (slope and intercept)
+/// \author Joshua König
+/// \ingroup EMCALCalib
+/// \since June 30, 2025
+
+#ifndef EMCALTEMPCALIBEXTRACTOR_H_
+#define EMCALTEMPCALIBEXTRACTOR_H_
+
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include "CCDB/BasicCCDBManager.h"
+#include "EMCALCalib/ElmbData.h"
+#include "EMCALCalib/TempCalibrationParams.h"
+#include "EMCALBase/Geometry.h"
+
+namespace o2
+{
+namespace emcal
+{
+
+class EMCALTempCalibExtractor
+{
+
+ public:
+  /// \brief Constructor
+  EMCALTempCalibExtractor()
+  {
+    LOG(info) << "initialized EMCALTempCalibExtractor";
+    try {
+      // Try to access geometry initialized ountside
+      mGeometry = o2::emcal::Geometry::GetInstance();
+    } catch (o2::emcal::GeometryNotInitializedException& e) {
+      mGeometry = o2::emcal::Geometry::GetInstanceFromRunNumber(300000); // fallback option
+    }
+  };
+  /// \brief Destructor
+  ~EMCALTempCalibExtractor() = default;
+
+  /// \brief Initialize temperature data and slope for each cell from the ccdb
+  /// \param path path to the slope data
+  /// \param timestamp timestamp for the ccdb objects or runnumber (will detect automatically if its a runnumber and convert it)
+  void InitializeFromCCDB(std::string path, uint64_t timestamp);
+
+  /// \brief get average temperature in a supermodule
+  /// \param iSM SM number
+  /// \param ElmbData object where temperature sensor values are stored
+  /// \return average temperature in a supermodule
+  float getTemperatureForSM(const unsigned short iSM, o2::emcal::ElmbData* ElmbData) const;
+
+  /// \brief get gain calibration factor depending on the temperature and the slope of the cell
+  /// \param cellID cell ID
+  /// \return gain calibration factor
+  float getGainCalibFactor(const unsigned short cellID) const;
+
+  /// \brief set temperature range in which sensor ddata is assumed to be good
+  /// \param low lower temperature
+  /// \param high upper temperature
+  void setAcceptedEnergyRange(float low, float high);
+
+  /// \brief set if median (true) or mean (false) should be used for averaging of the temperature in a SM
+  void setUseMedian(const bool tmp) { mUseMedian = tmp; }
+
+  /// \brief get sensor IDs for a specific supermodule
+  /// \param iSM SM number
+  /// \return vector of sensor IDs
+  std::vector<unsigned short> getSensorsForSM(const unsigned short iSM) const;
+
+ private:
+  static constexpr unsigned short mNCells = 17664;      ///< Number of EMCal cells
+  std::array<float, mNCells> mGainCalibFactors;         ///< gain calibration factors that are calculated based on the temperature and the slopes for each cell
+  o2::emcal::Geometry* mGeometry;                       ///< pointer to the EMCal geometry
+  std::array<float, 2> mAcceptedTempRange = {15., 30.}; ///< Temperature range where sensors are believed to send good data. Temperatures outside this range will be rejected
+  bool mUseMedian = true;                               /// switch to decide if temperature within a SM should be calculated as the mean or the median of the individual sensor data
+};
+
+} // namespace emcal
+
+} // namespace o2
+
+#endif