ITS3: harmonize study and alignment wfx for misalignment

Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>

Author: f3sch

Detectors/Upgrades/ITS3/alignment/CMakeLists.txt

@@ -13,6 +13,9 @@
 o2_add_library(ITS3Align
                TARGETVARNAME targetName
                SOURCES src/AlignmentHierarchy.cxx
+                       src/AlignmentDOF.cxx
+                       src/AlignmentMath.cxx
+                       src/MisalignmentUtils.cxx
                        src/AlignmentSensors.cxx
                        src/AlignmentParams.cxx
                        src/AlignmentTypes.cxx

Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentDOF.h

@@ -12,11 +12,29 @@
 #ifndef O2_ITS3_ALIGNMENT_DOF_H
 #define O2_ITS3_ALIGNMENT_DOF_H
 
+#include <algorithm>
 #include <cstdint>
 #include <format>
-#include <algorithm>
+#include <stdexcept>
+#include <string>
 #include <vector>
 
+#include <Eigen/Dense>
+
+struct DerivativeContext {
+  int sensorID{-1};
+  int layerID{-1};
+  double measX{0.};
+  double measAlpha{0.};
+  double measZ{0.};
+  double trkY{0.};
+  double trkZ{0.};
+  double snp{0.};
+  double tgl{0.};
+  double dydx{0.};
+  double dzdx{0.};
+};
+
 // Generic set of DOF
 class DOFSet
 {
@@ -30,6 +48,7 @@ class DOFSet
   virtual Type type() const = 0;
   int nDOFs() const { return static_cast<int>(mFree.size()); }
   virtual std::string dofName(int idx) const = 0;
+  virtual void fillDerivatives(const DerivativeContext& ctx, Eigen::Ref<Eigen::MatrixXd> out) const = 0;
   bool isFree(int idx) const { return mFree[idx]; }
   void setFree(int idx, bool f) { mFree[idx] = f; }
   void setAllFree(bool f) { std::fill(mFree.begin(), mFree.end(), f); }
@@ -73,6 +92,7 @@ class RigidBodyDOFSet final : public DOFSet
   }
   Type type() const override { return Type::RigidBody; }
   std::string dofName(int idx) const override { return RigidBodyDOFNames[idx]; }
+  void fillDerivatives(const DerivativeContext& ctx, Eigen::Ref<Eigen::MatrixXd> out) const override;
   uint8_t mask() const
   {
     uint8_t m = 0;
@@ -100,6 +120,7 @@ class LegendreDOFSet final : public DOFSet
     int j = idx - (i * (i + 1) / 2);
     return std::format("L({},{})", i, j);
   }
+  void fillDerivatives(const DerivativeContext& ctx, Eigen::Ref<Eigen::MatrixXd> out) const override;
 
  private:
   int mOrder;
@@ -146,6 +167,7 @@ class InextensionalDOFSet final : public DOFSet
     static constexpr const char* subNames[] = {"a", "b", "c", "d"};
     return std::format("{}_{}", subNames[sub], n);
   }
+  void fillDerivatives(const DerivativeContext& ctx, Eigen::Ref<Eigen::MatrixXd> out) const override;
 
  private:
   int mMaxOrder;

Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentHierarchy.h

@@ -33,24 +33,6 @@ namespace o2::its3::align
 using Matrix36 = Eigen::Matrix<double, 3, 6>;
 using Matrix66 = Eigen::Matrix<double, 6, 6>;
 
-// return the rigid body derivatives
-// trk has be at in the measurment frame
-auto getRigidBodyDerivatives(const auto& trk)
-{
-  // calculate slopes
-  const double tgl = trk.getTgl(), snp = trk.getSnp();
-  const double csp = 1. / sqrt(1. + (tgl * tgl));
-  const double u = trk.getY(), v = trk.getZ();
-  const double uP = snp * csp, vP = tgl * csp;
-  Matrix36 der;
-  der.setZero();
-  // columns: Tt,  Tu,  Tv,  Rt,    Ru,   Rv
-  //          (X)  (Y)  (Z)  (RX)   (RY)  (RZ)
-  der << uP, -1., 0., v, v * uP, -u * uP,
-    vP, 0., -1., -u, v * vP, -u * vP;
-  return der;
-}
-
 class HierarchyConstraint
 {
  public:

Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentMath.h

@@ -0,0 +1,32 @@
+// Copyright 2019-2026 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.
+
+#ifndef O2_ITS3_ALIGNMENT_MATH_H
+#define O2_ITS3_ALIGNMENT_MATH_H
+
+#include <utility>
+#include <vector>
+
+namespace o2::its3::align
+{
+
+struct TrackSlopes {
+  double dydx{0.};
+  double dzdx{0.};
+};
+
+std::pair<double, double> computeUV(double gloX, double gloY, double gloZ, int sensorID, double radius);
+TrackSlopes computeTrackSlopes(double snp, double tgl);
+std::vector<double> legendrePols(int order, double x);
+
+} // namespace o2::its3::align
+
+#endif

Detectors/Upgrades/ITS3/alignment/include/ITS3Align/MisalignmentUtils.h

@@ -0,0 +1,79 @@
+// Copyright 2019-2026 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.
+
+#ifndef O2_ITS3_ALIGNMENT_MISALIGNMENTUTILS_H
+#define O2_ITS3_ALIGNMENT_MISALIGNMENTUTILS_H
+
+#include <array>
+#include <cstddef>
+#include <map>
+#include <string>
+#include <utility>
+
+#include "ITS3Align/AlignmentMath.h"
+#include "MathUtils/LegendrePols.h"
+
+namespace o2::its3::align
+{
+
+struct InextensionalMisalignment {
+  std::map<int, std::array<double, 4>> modes; // n -> (a_n, b_n, c_n, d_n)
+  double alpha{0.};
+  double beta{0.};
+};
+
+struct SensorMisalignment {
+  o2::math_utils::Legendre2DPolynominal legendre;
+  bool hasLegendre{false};
+  InextensionalMisalignment inextensional;
+  bool hasInextensional{false};
+
+  bool empty() const noexcept { return !hasLegendre && !hasInextensional; }
+};
+
+struct MisalignmentModel {
+  static constexpr std::size_t NSensors = 6;
+  std::array<SensorMisalignment, NSensors> sensors{};
+
+  bool empty() const noexcept;
+  const SensorMisalignment& operator[](std::size_t idx) const { return sensors[idx]; }
+  SensorMisalignment& operator[](std::size_t idx) { return sensors[idx]; }
+};
+
+struct MisalignmentFrame {
+  int sensorID{-1};
+  int layerID{-1};
+  double x{0.};     // tracking-frame X / nominal radius at the measurement
+  double alpha{0.}; // tracking-frame alpha
+  double z{0.};     // tracking-frame measurement z
+};
+
+struct MisalignmentShift {
+  double dy{0.};
+  double dz{0.};
+  bool accepted{true};
+
+  MisalignmentShift& operator+=(const MisalignmentShift& other)
+  {
+    dy += other.dy;
+    dz += other.dz;
+    accepted = accepted && other.accepted;
+    return *this;
+  }
+};
+
+MisalignmentModel loadMisalignmentModel(const std::string& jsonPath);
+MisalignmentShift evaluateLegendreShift(const SensorMisalignment& sensor, const MisalignmentFrame& frame, const TrackSlopes& slopes);
+MisalignmentShift evaluateInextensionalShift(const SensorMisalignment& sensor, const MisalignmentFrame& frame, const TrackSlopes& slopes);
+
+} // namespace o2::its3::align
+
+#endif

Detectors/Upgrades/ITS3/alignment/src/AlignmentDOF.cxx

@@ -0,0 +1,114 @@
+// Copyright 2019-2026 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.
+
+#include "ITS3Align/AlignmentDOF.h"
+
+#include <cmath>
+#include <stdexcept>
+
+#include "ITS3Align/AlignmentMath.h"
+#include "ITS3Base/SpecsV2.h"
+
+namespace
+{
+
+void validateDerivativeOutput(const DOFSet& dofSet, Eigen::Ref<Eigen::MatrixXd> out)
+{
+  if (out.rows() != 3 || out.cols() != dofSet.nDOFs()) {
+    throw std::invalid_argument(std::format("Derivative buffer shape {}x{} does not match expected 3x{}",
+                                            out.rows(), out.cols(), dofSet.nDOFs()));
+  }
+  out.setZero();
+}
+
+} // namespace
+
+void RigidBodyDOFSet::fillDerivatives(const DerivativeContext& ctx, Eigen::Ref<Eigen::MatrixXd> out) const
+{
+  validateDerivativeOutput(*this, out);
+
+  const double csp = 1. / std::sqrt(1. + (ctx.tgl * ctx.tgl));
+  const double uP = ctx.snp * csp;
+  const double vP = ctx.tgl * csp;
+
+  out(0, TX) = uP;
+  out(0, TY) = -1.;
+  out(0, RX) = ctx.trkZ;
+  out(0, RY) = ctx.trkZ * uP;
+  out(0, RZ) = -ctx.trkY * uP;
+
+  out(1, TX) = vP;
+  out(1, TZ) = -1.;
+  out(1, RX) = -ctx.trkY;
+  out(1, RY) = ctx.trkZ * vP;
+  out(1, RZ) = -ctx.trkY * vP;
+}
+
+void LegendreDOFSet::fillDerivatives(const DerivativeContext& ctx, Eigen::Ref<Eigen::MatrixXd> out) const
+{
+  validateDerivativeOutput(*this, out);
+  if (ctx.sensorID < 0 || ctx.layerID < 0) {
+    throw std::invalid_argument("LegendreDOFSet requires an ITS3 measurement context");
+  }
+
+  const double gloX = ctx.measX * std::cos(ctx.measAlpha);
+  const double gloY = ctx.measX * std::sin(ctx.measAlpha);
+  const auto [u, v] = o2::its3::align::computeUV(gloX, gloY, ctx.measZ, ctx.sensorID, o2::its3::constants::radii[ctx.layerID]);
+  const auto pu = o2::its3::align::legendrePols(mOrder, u);
+  const auto pv = o2::its3::align::legendrePols(mOrder, v);
+
+  int idx = 0;
+  for (int i = 0; i <= mOrder; ++i) {
+    for (int j = 0; j <= i; ++j) {
+      const double basis = pu[j] * pv[i - j];
+      out(0, idx) = ctx.dydx * basis;
+      out(1, idx) = ctx.dzdx * basis;
+      ++idx;
+    }
+  }
+}
+
+void InextensionalDOFSet::fillDerivatives(const DerivativeContext& ctx, Eigen::Ref<Eigen::MatrixXd> out) const
+{
+  validateDerivativeOutput(*this, out);
+  if (ctx.layerID < 0) {
+    throw std::invalid_argument("InextensionalDOFSet requires an ITS3 measurement context");
+  }
+
+  const double r = o2::its3::constants::radii[ctx.layerID];
+  const double phi = std::atan2(r * std::sin(ctx.measAlpha), r * std::cos(ctx.measAlpha));
+  const double z = ctx.measZ;
+
+  for (int n = 2; n <= mMaxOrder; ++n) {
+    const double sn = std::sin(n * phi);
+    const double cn = std::cos(n * phi);
+    const double n2 = static_cast<double>(n * n);
+    const int off = modeOffset(n);
+
+    out(0, off + 0) = -(z / r) * (n * sn + ctx.dydx * n2 * cn);
+    out(1, off + 0) = -cn - ctx.dzdx * (z / r) * n2 * cn;
+
+    out(0, off + 1) = (z / r) * (n * cn - ctx.dydx * n2 * sn);
+    out(1, off + 1) = -sn * (1. + ctx.dzdx * (z / r) * n2);
+
+    out(0, off + 2) = -cn + ctx.dydx * n * sn;
+    out(1, off + 2) = ctx.dzdx * n * sn;
+
+    out(0, off + 3) = -sn - ctx.dydx * n * cn;
+    out(1, off + 3) = -ctx.dzdx * n * cn;
+  }
+
+  out(0, alphaIdx()) = z / r;
+  out(1, alphaIdx()) = -phi;
+
+  out(0, betaIdx()) = -phi - ctx.dydx;
+  out(1, betaIdx()) = -ctx.dzdx;
+}

Detectors/Upgrades/ITS3/alignment/src/AlignmentMath.cxx

@@ -0,0 +1,54 @@
+// Copyright 2019-2026 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.
+
+#include "ITS3Align/AlignmentMath.h"
+
+#include <cmath>
+
+#include <TMath.h>
+
+#include "ITS3Base/SpecsV2.h"
+#include "MathUtils/Utils.h"
+
+namespace o2::its3::align
+{
+
+std::pair<double, double> computeUV(double gloX, double gloY, double gloZ, int sensorID, double radius)
+{
+  const bool isTop = sensorID % 2 == 0;
+  const double phi = o2::math_utils::to02Pid(std::atan2(gloY, gloX));
+  const double phiBorder1 = o2::math_utils::to02Pid(((isTop ? 0. : 1.) * TMath::Pi()) + std::asin(constants::equatorialGap / 2. / radius));
+  const double phiBorder2 = o2::math_utils::to02Pid(((isTop ? 1. : 2.) * TMath::Pi()) - std::asin(constants::equatorialGap / 2. / radius));
+  const double u = (((phi - phiBorder1) * 2.) / (phiBorder2 - phiBorder1)) - 1.;
+  const double v = ((2. * gloZ + constants::segment::lengthSensitive) / constants::segment::lengthSensitive) - 1.;
+  return {u, v};
+}
+
+TrackSlopes computeTrackSlopes(double snp, double tgl)
+{
+  const double csci = 1. / std::sqrt(1. - (snp * snp));
+  return {.dydx = snp * csci, .dzdx = tgl * csci};
+}
+
+std::vector<double> legendrePols(int order, double x)
+{
+  std::vector<double> p(order + 1);
+  p[0] = 1.;
+  if (order > 0) {
+    p[1] = x;
+  }
+  for (int n = 1; n < order; ++n) {
+    p[n + 1] = ((2 * n + 1) * x * p[n] - n * p[n - 1]) / (n + 1);
+  }
+  return p;
+}
+
+} // namespace o2::its3::align