ITS3: fix inextensional model + opt. radial

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

Author: f3sch

Detectors/Upgrades/ITS3/alignment/README.md

@@ -31,7 +31,23 @@ dofSet.json:
 ```
 
 
-## In-existensional modes
+## In-extensional modes
+
+The deformation of the open half-shell is parameterised by two 1D functions expanded in Legendre polynomials of the
+normalised azimuth `u` (the same coordinate as the radial Legendre model):
+
+```
+f(phi) = sum_k f_k P_k(u),   g(phi) = sum_k g_k P_k(u)
+u_z = f,   u_phi = -(z/r) f' + g,   u_r = (z/r) f'' - g'
+```
+
+`order` sets the maximum `k`. Optionally, strictly radial ("extensional") modes can be added on top, `u_r += sum_{k,l}
+h_{k,l} P_k(u) P_l(v)` with `l >= 1`, enabled via `extOrderPhi` (max `k`) and `extOrderZ` (max `l`); `l = 0` is excluded
+because a z-independent radial field is already spanned by the `g` family.
+
+Note that `f_0` (translation along the cylinder axis) and `g_0` (rotation about it) are rigid-body motions and are fixed
+by default; free them only if the rigid-body DOFs of the same volume are not fitted.
+
 ```json
 {
   "defaults": { "rigidBody": "fixed" },
@@ -40,24 +56,26 @@ dofSet.json:
       "match": "ITS3Layer1/ITS3CarbonForm0",
       "calib": {
         "type": "inextensional",
-        "order": 2,
-        "free": ["a_2", "b_2", "c_2", "d_2", "alpha", "beta"]
+        "order": 10,
+        "extOrderPhi": 7,
+        "extOrderZ": 8,
+        "fix": ["f_0", "g_0"]
       }
     }
   ]
 }
 ```
 
+Injected/fitted coefficients (`h` keys are `"<k>_<l>"`):
+
 ```json
 [
   {
     "id": 2,
     "inextensional": {
-      "modes": {
-        "2": [0.0008, -0.0005, 0.0006, -0.0007]
-      },
-      "alpha": 0.0004,
-      "beta": -0.0003
+      "f": { "1": 0.0001, "2": -0.0002 },
+      "g": { "1": 0.0625, "3": 0.0335, "5": -0.0453 },
+      "h": { "4_2": -0.0421, "6_2": 0.0252, "4_4": 0.0435 }
     }
   }
 ]

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

@@ -21,6 +21,8 @@
 
 #include <Eigen/Dense>
 
+#include "ITS3Align/AlignmentLabel.h"
+
 struct DerivativeContext {
   int sensorID{-1};
   int layerID{-1};
@@ -62,7 +64,7 @@ class DOFSet
   }
 
  protected:
-  DOFSet(int n) : mFree(n, true) {}
+  DOFSet(int n) : mFree(n, true) { GlobalLabel::checkDOFCount(n); }
   std::vector<bool> mFree;
 };
 
@@ -126,51 +128,77 @@ class LegendreDOFSet final : public DOFSet
   int mOrder;
 };
 
-// In-extensional deformation DOFs for cylindrical half-shells
-// Fourier modes n=2..N: 4 params each (a_n, b_n, c_n, d_n)
-// Plus 2 non-periodic modes (alpha, beta) for the half-cylinder open edges
-// Total: 4*(N-1) + 2
+// Deformation DOFs for an open cylindrical half-shell.
+//
+// Inextensional part. Vanishing linear membrane strains admit the general solution (u in the local (r, phi, z)
+// directions) u_z = f(phi) u_phi = -(z/r) f'(phi) + g(phi) u_r = (z/r) f''(phi) - g'(phi) with two arbitrary
+// one-dimensional functions f, g. Because the shell is open in phi these are expanded in Legendre polynomials of the
+// normalised azimuth u in [-1, 1]: f(phi) = sum_k f_k P_k(u), g(phi) = sum_k g_k P_k(u).
+//
+// Extensional part (optional). The inextensional u_r is at most linear in z, so radial deformations with curvature
+// along z lie outside it. They are added as strictly radial modes u_r += sum_{k,l} h_{k,l} P_k(u) P_l(v), l >= 1, with
+// v the normalised axial coordinate. l = 0 is excluded because a z-independent radial field is already spanned by the g
+// family.
+//
+// Flat index layout: [f_0, g_0, f_1, g_1, ..., f_K, g_K, h_{0,1} ... h_{0,Lz}, h_{1,1} ... h_{Kphi,Lz}]
+//
+// NOTE on degeneracies: f_0 is a rigid translation along the cylinder axis and g_0 a rigid rotation about it, i.e. they
+// duplicate rigid-body DOFs of the same volume.
 class InextensionalDOFSet final : public DOFSet
 {
  public:
-  explicit InextensionalDOFSet(int maxOrder) : DOFSet((4 * (maxOrder - 1)) + 2), mMaxOrder(maxOrder)
+  explicit InextensionalDOFSet(int maxOrder, int extOrderPhi = -1, int extOrderZ = 0)
+    : DOFSet(nDOFsFor(maxOrder, extOrderPhi, extOrderZ)),
+      mMaxOrder(maxOrder),
+      mExtOrderPhi(extOrderZ > 0 ? extOrderPhi : -1),
+      mExtOrderZ(extOrderPhi >= 0 ? extOrderZ : 0)
   {
-    if (maxOrder < 2) {
-      // the rest is eq. to rigid body
-      throw std::invalid_argument("InextensionalDOFSet requires maxOrder >= 2");
+    if (maxOrder < 1) {
+      // only k = 0 is left, which is equivalent to a rigid body motion
+      throw std::invalid_argument("InextensionalDOFSet requires maxOrder >= 1");
     }
+    // f_0 / g_0 are rigid: fixed unless explicitly freed
+    setFree(fIdx(0), false);
+    setFree(gIdx(0), false);
   }
+
+  static int nDOFsFor(int maxOrder, int extOrderPhi, int extOrderZ)
+  {
+    int n = 2 * (maxOrder + 1);
+    if (extOrderPhi >= 0 && extOrderZ > 0) {
+      n += (extOrderPhi + 1) * extOrderZ;
+    }
+    return n;
+  }
+
   Type type() const override { return Type::Inextensional; }
   int maxOrder() const { return mMaxOrder; }
+  int extOrderPhi() const { return mExtOrderPhi; }
+  int extOrderZ() const { return mExtOrderZ; }
+  bool hasExtensional() const { return mExtOrderPhi >= 0 && mExtOrderZ > 0; }
 
-  // number of periodic DOFs (before alpha, beta)
-  int nPeriodic() const { return 4 * (mMaxOrder - 1); }
-
-  // flat index layout: [a_2, b_2, c_2, d_2, a_3, b_3, c_3, d_3, ..., alpha, beta]
-  // index of first DOF for mode n
-  static int modeOffset(int n) { return 4 * (n - 2); }
+  // number of inextensional DOFs (before the radial h modes)
+  int nInextensional() const { return 2 * (mMaxOrder + 1); }
 
-  // indices of the non-periodic modes
-  int alphaIdx() const { return nPeriodic(); }
-  int betaIdx() const { return nPeriodic() + 1; }
+  // flat indices
+  static int fIdx(int k) { return 2 * k; }
+  static int gIdx(int k) { return (2 * k) + 1; }
+  int hIdx(int k, int l) const { return nInextensional() + (k * mExtOrderZ) + (l - 1); }
 
   std::string dofName(int idx) const override
   {
-    if (idx == alphaIdx()) {
-      return "alpha";
-    }
-    if (idx == betaIdx()) {
-      return "beta";
+    if (idx < nInextensional()) {
+      return std::format("{}_{}", (idx % 2 == 0) ? "f" : "g", idx / 2);
     }
-    int n = (idx / 4) + 2;
-    int sub = idx % 4;
-    static constexpr const char* subNames[] = {"a", "b", "c", "d"};
-    return std::format("{}_{}", subNames[sub], n);
+    const int e = idx - nInextensional();
+    return std::format("h_{}_{}", e / mExtOrderZ, (e % mExtOrderZ) + 1);
   }
   void fillDerivatives(const DerivativeContext& ctx, Eigen::Ref<Eigen::MatrixXd> out) const override;
 
  private:
   int mMaxOrder;
+  int mExtOrderPhi;
+  int mExtOrderZ;
 };
 
 #endif

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

@@ -13,18 +13,19 @@
 #define O2_ITS3_ALIGNMENT_LABEL_H
 
 #include <cstdint>
+#include <stdexcept>
 #include <string>
 #include <format>
 
 class GlobalLabel
 {
   // Millepede label is any positive integer [1....)
-  // Layout: DOF(5) | CALIB(1) | ID(22) | SENS(1) | DET(2) = 31 usable bits (MSB reserved, GBL uses signed int)
+  // Layout: DOF(8) | CALIB(1) | ID(19) | SENS(1) | DET(2) = 31 usable bits (MSB reserved, GBL uses signed int)
  public:
   using T = uint32_t;
-  static constexpr int DOF_BITS = 5;   // bits 0-4
-  static constexpr int CALIB_BITS = 1; // bit 5: 0 = rigid body, 1 = calibration (only allow for one calibration, could be extended if needed)
-  static constexpr int ID_BITS = 22;   // bits 6-27
+  static constexpr int DOF_BITS = 8;   // bits 0-7
+  static constexpr int CALIB_BITS = 1; // bit 8: 0 = rigid body, 1 = calibration (only allow for one calibration, could be extended if needed)
+  static constexpr int ID_BITS = 19;   // bits 9-27
   static constexpr int SENS_BITS = 1;  // bit 28
   static constexpr int TOTAL_BITS = sizeof(T) * 8;
   static constexpr int DET_BITS = TOTAL_BITS - (DOF_BITS + CALIB_BITS + ID_BITS + SENS_BITS) - 1; // one less bit since GBL uses int!
@@ -48,12 +49,30 @@ class GlobalLabel
   static constexpr T DET_MAX = (T(1) << DET_BITS) - T(1);
   static constexpr T DET_MASK = DET_MAX << DET_SHIFT;
 
+  /// maximum number of DOFs that can be labelled on one volume (per calib bit)
+  static constexpr int MAX_DOFS = static_cast<int>(DOF_MAX) + 1;
+
+  /// throws if a DOF set is too large to be labelled without aliasing
+  static void checkDOFCount(int nDOFs)
+  {
+    if (nDOFs > MAX_DOFS) {
+      throw std::out_of_range(std::format(
+        "DOF set with {} parameters exceeds the {} labelable DOFs (DOF_BITS={}); "
+        "distinct parameters would alias onto the same Millepede label",
+        nDOFs, MAX_DOFS, DOF_BITS));
+    }
+  }
+
   GlobalLabel(T det, T id, bool sens, bool calib = false)
     : mID((((id + 1) & ID_MAX) << ID_SHIFT) |
           ((det & DET_MAX) << DET_SHIFT) |
           ((T(sens) & SENS_MAX) << SENS_SHIFT) |
           ((T(calib) & CALIB_MAX) << CALIB_SHIFT))
   {
+    if ((id + 1) > ID_MAX) {
+      throw std::out_of_range(std::format("Volume id {} exceeds the {} labelable ids (ID_BITS={})",
+                                          id, ID_MAX - 1, ID_BITS));
+    }
   }
 
   /// produce the raw Millepede label for a given DOF index (rigid body: calib=0 in label)

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

@@ -23,11 +23,18 @@ struct TrackSlopes {
   double dzdx{0.};
 };
 
-double getSensorPhiWidth(int sensorID, double radius);
 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);
 
+// First and second derivatives dP_n/dx, d^2P_n/dx^2 for n = 0..order.
+std::vector<double> legendrePolsD1(int order, double x);
+std::vector<double> legendrePolsD2(int order, double x);
+
+// Jacobian factor of the angular normalisation used by computeUV: c_phi = du/dphi = 2 / (phiBorder2 - phiBorder1).
+// Needed to convert derivatives with respect to the normalised u back to derivatives with respect to the azimuth phi.
+double phiScale(double radius);
+
 } // namespace o2::its3::align
 
 #endif

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

@@ -24,10 +24,13 @@
 namespace o2::its3::align
 {
 
+// Legendre parameterisation of the open half-shell deformation, matching InextensionalDOFSet: inextensional
+// coefficients f_k, g_k of the normalised azimuth u, plus optional strictly radial modes h_{k,l} (l >= 1) in P_k(u)
+// P_l(v). See AlignmentDOF.h for the displacement field.
 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.};
+  std::map<int, double> f;                 // k -> f_k
+  std::map<int, double> g;                 // k -> g_k
+  std::map<std::pair<int, int>, double> h; // (k, l) -> h_{k,l}, l >= 1
 };
 
 struct SensorMisalignment {