ITS3: alignment account for tangential movement

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

Author: f3sch

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

@@ -23,6 +23,7 @@ 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);

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

@@ -11,11 +11,13 @@
 
 #include "ITS3Align/AlignmentDOF.h"
 
+#include <array>
 #include <cmath>
 #include <stdexcept>
 
 #include "ITS3Align/AlignmentMath.h"
 #include "ITS3Base/SpecsV2.h"
+#include "CommonConstants/MathConstants.h"
 
 namespace
 {
@@ -64,12 +66,34 @@ void LegendreDOFSet::fillDerivatives(const DerivativeContext& ctx, Eigen::Ref<Ei
   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);
+  const double phiWidth = o2::its3::align::getSensorPhiWidth(ctx.sensorID, o2::its3::constants::radii[ctx.layerID]);
+
+  // same intergration as `evaluateLegendreShift' but now for each order separateley
+  Eigen::VectorXd arcMismatch = Eigen::VectorXd::Zero(nDOFs());
+  if (std::abs(u) > o2::constants::math::Almost0) {
+    constexpr std::array<double, 8> x = {-0.9602898564975363, -0.7966664774136267, -0.5255324099163290, -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, 0.7966664774136267, 0.9602898564975363};
+    constexpr std::array<double, 8> w = {0.1012285362903763, 0.2223810344533745, 0.3137066458778873, 0.3626837833783620, 0.3626837833783620, 0.3137066458778873, 0.2223810344533745, 0.1012285362903763};
+    const double mid = 0.5 * u;
+    const double half = 0.5 * u;
+    for (int iq = 0; iq < 8; ++iq) {
+      const double up = mid + (half * x[iq]);
+      const auto puQ = o2::its3::align::legendrePols(mOrder, up);
+      int idx = 0;
+      for (int i = 0; i <= mOrder; ++i) {
+        for (int j = 0; j <= i; ++j) {
+          arcMismatch[idx] += w[iq] * puQ[j] * pv[i - j];
+          ++idx;
+        }
+      }
+    }
+    arcMismatch *= 0.5 * phiWidth * half;
+  }
 
   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(0, idx) = (ctx.dydx * basis) + arcMismatch[idx];
       out(1, idx) = ctx.dzdx * basis;
       ++idx;
     }

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

@@ -21,6 +21,15 @@
 namespace o2::its3::align
 {
 
+double getSensorPhiWidth(int sensorID, double radius)
+{
+  const bool isTop = sensorID % 2 == 0;
+  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 width = phiBorder2 - phiBorder1;
+  return (width < 0.) ? width + TMath::TwoPi() : width;
+}
+
 std::pair<double, double> computeUV(double gloX, double gloY, double gloZ, int sensorID, double radius)
 {
   const bool isTop = sensorID % 2 == 0;

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

@@ -206,6 +206,10 @@ void AlignmentSpec::process()
     buildT2V();
   }
 
+  if (mNThreads > 1 && !(mParams->misAlgJson.empty())) {
+    LOGP(warn, "Applying misalignment works only single-threaded, forcing to 1");
+    mNThreads = 1;
+  }
   LOGP(info, "Starting fits with {} threads", mNThreads);
 
   // Data
@@ -831,6 +835,9 @@ bool AlignmentSpec::applyMisalignment(Eigen::Vector2d& res, const FrameInfoExt&
     const auto prop = o2::base::PropagatorD::Instance();
 
     const auto lbl = mRecoData->getITSTracksMCLabels()[iTrk];
+    if (lbl.isFake()) {
+      return false;
+    }
     const auto mcTrk = mcReader->getTrack(lbl);
     if (!mcTrk) {
       return false;

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

@@ -15,11 +15,13 @@
 #include <cmath>
 #include <fstream>
 #include <vector>
+#include <array>
 
 #include <TMatrixD.h>
 #include <nlohmann/json.hpp>
 
 #include "Framework/Logger.h"
+#include "CommonConstants/MathConstants.h"
 #include "ITS3Base/SpecsV2.h"
 
 namespace o2::its3::align
@@ -98,9 +100,27 @@ MisalignmentShift evaluateLegendreShift(const SensorMisalignment& sensor, const
   auto [u, v] = computeUV(gloX, gloY, gloZ, frame.sensorID, constants::radii[frame.layerID]);
   const double h = sensor.legendre(u, v);
 
+  // this is the shift due to back-projection of the track on the ideal surface
   shift.dy = slopes.dydx * h;
   shift.dz = slopes.dzdx * h;
 
+  if (std::abs(u) > o2::constants::math::Almost0) {
+    // account for additional tangential movement due to radial shift
+    // we have to approximate the difference in arc-length from the reference pnt on the deformed surface
+    // this is done by integrating the height function via Gauss-Legendre quadrature (from Numerical recipes 4.6 [1])
+    constexpr std::array<double, 8> x = {-0.9602898564975363, -0.7966664774136267, -0.5255324099163290, -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, 0.7966664774136267, 0.9602898564975363};
+    constexpr std::array<double, 8> w = {0.1012285362903763, 0.2223810344533745, 0.3137066458778873, 0.3626837833783620, 0.3626837833783620, 0.3137066458778873, 0.2223810344533745, 0.1012285362903763};
+    const double mid = 0.5 * u;
+    const double half = 0.5 * u;
+    double integral = 0.;
+    for (int i = 0; i < 8; ++i) {
+      const double up = mid + (half * x[i]);
+      integral += w[i] * sensor.legendre(up, v);
+    }
+    integral *= half;
+    shift.dy += 0.5 * getSensorPhiWidth(frame.sensorID, constants::radii[frame.layerID]) * integral;
+  }
+
   const double newGloY = gloY + (shift.dy * std::cos(frame.alpha));
   const double newGloX = gloX - (shift.dy * std::sin(frame.alpha));
   const double newGloZ = gloZ + shift.dz;