From 1a4325ec2aa76ad47ac107bace6842529b6b4c88 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 17 Jun 2026 11:24:35 +0200 Subject: [PATCH] ITS3: fix inextensional model + opt. radial Signed-off-by: Felix Schlepper --- Detectors/Upgrades/ITS3/alignment/README.md | 34 +++-- .../include/ITS3Align/AlignmentDOF.h | 82 +++++++---- .../include/ITS3Align/AlignmentLabel.h | 27 +++- .../include/ITS3Align/AlignmentMath.h | 9 +- .../include/ITS3Align/MisalignmentUtils.h | 9 +- .../ITS3/alignment/src/AlignmentDOF.cxx | 107 +++++++-------- .../ITS3/alignment/src/AlignmentHierarchy.cxx | 85 +++++++----- .../ITS3/alignment/src/AlignmentMath.cxx | 42 ++++-- .../ITS3/alignment/src/MisalignmentUtils.cxx | 116 +++++++++------- .../ITS3/study/macros/PlotMisalignment.C | 37 ++++- .../Upgrades/ITS3/study/src/TrackingStudy.cxx | 127 ++++++++++++------ 11 files changed, 429 insertions(+), 246 deletions(-) diff --git a/Detectors/Upgrades/ITS3/alignment/README.md b/Detectors/Upgrades/ITS3/alignment/README.md index 80213eb4e03b1..7d51df0465e44 100644 --- a/Detectors/Upgrades/ITS3/alignment/README.md +++ b/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 `"_"`): + ```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 } } } ] diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentDOF.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentDOF.h index 3fed9decbd6e7..991f074ed398f 100644 --- a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentDOF.h +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentDOF.h @@ -21,6 +21,8 @@ #include +#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 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 out) const override; private: int mMaxOrder; + int mExtOrderPhi; + int mExtOrderZ; }; #endif diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentLabel.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentLabel.h index 83495491b87e0..32ff02231ef78 100644 --- a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentLabel.h +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentLabel.h @@ -13,18 +13,19 @@ #define O2_ITS3_ALIGNMENT_LABEL_H #include +#include #include #include 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(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) diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentMath.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentMath.h index b3b538ab10192..c3aaa30cb5014 100644 --- a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentMath.h +++ b/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 computeUV(double gloX, double gloY, double gloZ, int sensorID, double radius); TrackSlopes computeTrackSlopes(double snp, double tgl); std::vector legendrePols(int order, double x); +// First and second derivatives dP_n/dx, d^2P_n/dx^2 for n = 0..order. +std::vector legendrePolsD1(int order, double x); +std::vector 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 diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/MisalignmentUtils.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/MisalignmentUtils.h index 457eccaeff4e6..e821161440ba1 100644 --- a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/MisalignmentUtils.h +++ b/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> modes; // n -> (a_n, b_n, c_n, d_n) - double alpha{0.}; - double beta{0.}; + std::map f; // k -> f_k + std::map g; // k -> g_k + std::map, double> h; // (k, l) -> h_{k,l}, l >= 1 }; struct SensorMisalignment { diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentDOF.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentDOF.cxx index 351cc0027f4c9..891f21fae5edf 100644 --- a/Detectors/Upgrades/ITS3/alignment/src/AlignmentDOF.cxx +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentDOF.cxx @@ -11,13 +11,11 @@ #include "ITS3Align/AlignmentDOF.h" -#include #include #include #include "ITS3Align/AlignmentMath.h" #include "ITS3Base/SpecsV2.h" -#include "CommonConstants/MathConstants.h" namespace { @@ -36,22 +34,17 @@ void validateDerivativeOutput(const DOFSet& dofSet, Eigen::Ref void RigidBodyDOFSet::fillDerivatives(const DerivativeContext& ctx, Eigen::Ref 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, TX) = ctx.dydx; out(0, TY) = -1.; out(0, RX) = ctx.trkZ; - out(0, RY) = ctx.trkZ * uP; - out(0, RZ) = -ctx.trkY * uP; + out(0, RY) = ctx.trkZ * ctx.dydx; + out(0, RZ) = -ctx.trkY * ctx.dydx; - out(1, TX) = vP; + out(1, TX) = ctx.dzdx; out(1, TZ) = -1.; out(1, RX) = -ctx.trkY; - out(1, RY) = ctx.trkZ * vP; - out(1, RZ) = -ctx.trkY * vP; + out(1, RY) = ctx.trkZ * ctx.dzdx; + out(1, RZ) = -ctx.trkY * ctx.dzdx; } void LegendreDOFSet::fillDerivatives(const DerivativeContext& ctx, Eigen::Ref out) const @@ -66,34 +59,12 @@ void LegendreDOFSet::fillDerivatives(const DerivativeContext& ctx, Eigen::Ref o2::constants::math::Almost0) { - constexpr std::array x = {-0.9602898564975363, -0.7966664774136267, -0.5255324099163290, -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, 0.7966664774136267, 0.9602898564975363}; - constexpr std::array 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) + arcMismatch[idx]; + out(0, idx) = ctx.dydx * basis; out(1, idx) = ctx.dzdx * basis; ++idx; } @@ -103,36 +74,48 @@ void LegendreDOFSet::fillDerivatives(const DerivativeContext& ctx, Eigen::Ref out) const { validateDerivativeOutput(*this, out); - if (ctx.layerID < 0) { + if (ctx.sensorID < 0 || 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(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; + 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, r); + const double cPhi = o2::its3::align::phiScale(r); + const double zOverR = ctx.measZ / r; + + // The residual derivative for a mode with displacement M = (M_r, M_phi, M_z) + // along the local (r, phi, z) directions is + // row0 = dydx * M_r - M_phi, row1 = dzdx * M_r - M_z + // (cf. the rigid-body case, where M = (1,0,0) gives (dydx, dzdx)). + const auto fill = [&out, &ctx](int idx, double mR, double mPhi, double mZ) { + out(0, idx) = (ctx.dydx * mR) - mPhi; + out(1, idx) = (ctx.dzdx * mR) - mZ; + }; + + const int order = std::max(mMaxOrder, hasExtensional() ? mExtOrderPhi : 0); + const auto pu = o2::its3::align::legendrePols(order, u); + const auto pu1 = o2::its3::align::legendrePolsD1(mMaxOrder, u); + const auto pu2 = o2::its3::align::legendrePolsD2(mMaxOrder, u); + + for (int k = 0; k <= mMaxOrder; ++k) { + // f_k: u_z = P_k, u_phi = -(z/r) c P'_k, u_r = (z/r) c^2 P''_k + fill(fIdx(k), + zOverR * cPhi * cPhi * pu2[k], + -zOverR * cPhi * pu1[k], + pu[k]); + // g_k: u_phi = P_k, u_r = -c P'_k + fill(gIdx(k), -cPhi * pu1[k], pu[k], 0.); } - out(0, alphaIdx()) = z / r; - out(1, alphaIdx()) = -phi; - - out(0, betaIdx()) = -phi - ctx.dydx; - out(1, betaIdx()) = -ctx.dzdx; + if (hasExtensional()) { + const auto pv = o2::its3::align::legendrePols(mExtOrderZ, v); + for (int k = 0; k <= mExtOrderPhi; ++k) { + for (int l = 1; l <= mExtOrderZ; ++l) { + // h_{k,l}: strictly radial u_r = P_k(u) P_l(v) + fill(hIdx(k, l), pu[k] * pv[l], 0., 0.); + } + } + } } diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx index 938c14c2c4759..8608248e60da8 100644 --- a/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx @@ -364,7 +364,10 @@ void applyDOFConfig(AlignableVolume* root, const std::string& jsonPath) vol->setCalib(std::move(dofSet)); } else if (calType == "inextensional") { int maxOrder = cal.value("order", 2); - auto dofSet = std::make_unique(maxOrder); + // optional strictly radial (extensional) modes h_{k,l}, l >= 1 + int extOrderPhi = cal.value("extOrderPhi", -1); + int extOrderZ = cal.value("extOrderZ", 0); + auto dofSet = std::make_unique(maxOrder, extOrderPhi, extOrderZ); bool fixed = cal.value("fixed", false); if (fixed) { dofSet->setAllFree(false); @@ -439,9 +442,9 @@ void writeMillepedeResults(AlignableVolume* root, const std::string& milleResPat std::map> injRB; std::map>> injMatrix; struct InjInex { - std::map> modes; - double alpha{0.}; - double beta{0.}; + std::map f; + std::map g; + std::map, double> h; }; std::map injInex; if (!injectedJsonPath.empty()) { @@ -459,16 +462,24 @@ void writeMillepedeResults(AlignableVolume* root, const std::string& milleResPat if (item.contains("inextensional")) { InjInex ii; const auto& inex = item["inextensional"]; - if (inex.contains("modes")) { - for (auto& [key, val] : inex["modes"].items()) { - ii.modes[std::stoi(key)] = val.get>(); + if (inex.contains("f")) { + for (auto& [key, val] : inex["f"].items()) { + ii.f[std::stoi(key)] = val.get(); } } - if (inex.contains("alpha")) { - ii.alpha = inex["alpha"].get(); + if (inex.contains("g")) { + for (auto& [key, val] : inex["g"].items()) { + ii.g[std::stoi(key)] = val.get(); + } } - if (inex.contains("beta")) { - ii.beta = inex["beta"].get(); + if (inex.contains("h")) { + for (auto& [key, val] : inex["h"].items()) { + const auto sep = key.find('_'); + if (sep == std::string::npos) { + continue; + } + ii.h[{std::stoi(key.substr(0, sep)), std::stoi(key.substr(sep + 1))}] = val.get(); + } } injInex[id] = ii; } @@ -533,38 +544,38 @@ void writeMillepedeResults(AlignableVolume* root, const std::string& milleResPat } else if (cal && cal->nFreeDOFs() && cal->type() == DOFSet::Type::Inextensional) { write = true; auto* inexSet = static_cast(cal); - int maxN = inexSet->maxOrder(); auto calibLbl = vol->getLabel().asCalib(); const auto& inj = injInex.contains(id) ? injInex[id] : InjInex{}; + auto fittedAt = [&](int idx) { + uint32_t raw = calibLbl.raw(idx); + auto it = labelToValue.find(raw); + return it != labelToValue.end() ? it->second : 0.0; + }; + json inexEntry; - json modesObj = json::object(); - for (int n = 2; n <= maxN; ++n) { - int off = InextensionalDOFSet::modeOffset(n); - std::array injCoeffs = {0., 0., 0., 0.}; - if (inj.modes.contains(n)) { - injCoeffs = inj.modes.at(n); - } - json modeArr = json::array(); - for (int k = 0; k < 4; ++k) { - uint32_t raw = calibLbl.raw(off + k); - auto it = labelToValue.find(raw); - double fitted = it != labelToValue.end() ? it->second : 0.0; - modeArr.push_back(fitted - injCoeffs[k]); - } - modesObj[std::to_string(n)] = modeArr; + json fObj = json::object(); + json gObj = json::object(); + for (int k = 0; k <= inexSet->maxOrder(); ++k) { + const double injF = inj.f.contains(k) ? inj.f.at(k) : 0.0; + const double injG = inj.g.contains(k) ? inj.g.at(k) : 0.0; + fObj[std::to_string(k)] = fittedAt(InextensionalDOFSet::fIdx(k)) - injF; + gObj[std::to_string(k)] = fittedAt(InextensionalDOFSet::gIdx(k)) - injG; } - inexEntry["modes"] = modesObj; - - // alpha - uint32_t rawAlpha = calibLbl.raw(inexSet->alphaIdx()); - auto itA = labelToValue.find(rawAlpha); - inexEntry["alpha"] = (itA != labelToValue.end() ? itA->second : 0.0) - inj.alpha; + inexEntry["f"] = fObj; + inexEntry["g"] = gObj; - // beta - uint32_t rawBeta = calibLbl.raw(inexSet->betaIdx()); - auto itB = labelToValue.find(rawBeta); - inexEntry["beta"] = (itB != labelToValue.end() ? itB->second : 0.0) - inj.beta; + if (inexSet->hasExtensional()) { + json hObj = json::object(); + for (int k = 0; k <= inexSet->extOrderPhi(); ++k) { + for (int l = 1; l <= inexSet->extOrderZ(); ++l) { + const auto key = std::pair{k, l}; + const double injH = inj.h.contains(key) ? inj.h.at(key) : 0.0; + hObj[std::format("{}_{}", k, l)] = fittedAt(inexSet->hIdx(k, l)) - injH; + } + } + inexEntry["h"] = hObj; + } entry["inextensional"] = inexEntry; } diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentMath.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentMath.cxx index 75dad6fb8af3d..e7093c394d0b4 100644 --- a/Detectors/Upgrades/ITS3/alignment/src/AlignmentMath.cxx +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentMath.cxx @@ -21,15 +21,6 @@ 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 computeUV(double gloX, double gloY, double gloZ, int sensorID, double radius) { const bool isTop = sensorID % 2 == 0; @@ -60,4 +51,37 @@ std::vector legendrePols(int order, double x) return p; } +std::vector legendrePolsD1(int order, double x) +{ + // P'_0 = 0, P'_1 = 1, P'_{n+1} = (2n+1) P_n + P'_{n-1} + const auto p = legendrePols(order, x); + std::vector d(order + 1, 0.); + if (order > 0) { + d[1] = 1.; + } + for (int n = 1; n < order; ++n) { + d[n + 1] = ((2 * n + 1) * p[n]) + d[n - 1]; + } + return d; +} + +std::vector legendrePolsD2(int order, double x) +{ + // P''_0 = P''_1 = 0, P''_{n+1} = (2n+1) P'_n + P''_{n-1} + const auto d1 = legendrePolsD1(order, x); + std::vector d(order + 1, 0.); + for (int n = 1; n < order; ++n) { + d[n + 1] = ((2 * n + 1) * d1[n]) + d[n - 1]; + } + return d; +} + +double phiScale(double radius) +{ + // computeUV maps [phiBorder1, phiBorder2] -> [-1, 1]; the span is the same + // for the top and bottom half-shells: pi - 2 asin(gap / 2r). + const double span = TMath::Pi() - (2. * std::asin(constants::equatorialGap / 2. / radius)); + return 2. / span; +} + } // namespace o2::its3::align diff --git a/Detectors/Upgrades/ITS3/alignment/src/MisalignmentUtils.cxx b/Detectors/Upgrades/ITS3/alignment/src/MisalignmentUtils.cxx index cd09c74c69ec9..cf5d0af417e5b 100644 --- a/Detectors/Upgrades/ITS3/alignment/src/MisalignmentUtils.cxx +++ b/Detectors/Upgrades/ITS3/alignment/src/MisalignmentUtils.cxx @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -21,7 +22,6 @@ #include #include "Framework/Logger.h" -#include "CommonConstants/MathConstants.h" #include "ITS3Base/SpecsV2.h" namespace o2::its3::align @@ -70,16 +70,41 @@ MisalignmentModel loadMisalignmentModel(const std::string& jsonPath) if (item.contains("inextensional")) { const auto& inex = item["inextensional"]; sensor.hasInextensional = true; - if (inex.contains("modes")) { - for (const auto& [key, val] : inex["modes"].items()) { - sensor.inextensional.modes[std::stoi(key)] = val.get>(); + // {"f": {"1": ..., "2": ...}, "g": {...}, "h": {"2_1": ..., "4_2": ...}} + if (inex.contains("f")) { + for (const auto& [key, val] : inex["f"].items()) { + sensor.inextensional.f[std::stoi(key)] = val.get(); } } - if (inex.contains("alpha")) { - sensor.inextensional.alpha = inex["alpha"].get(); + if (inex.contains("g")) { + for (const auto& [key, val] : inex["g"].items()) { + sensor.inextensional.g[std::stoi(key)] = val.get(); + } + } + if (inex.contains("h")) { + for (const auto& [key, val] : inex["h"].items()) { + const auto sep = key.find('_'); + if (sep == std::string::npos) { + LOGP(fatal, "Inextensional h key '{}' for sensor {} must be of the form '_' in {}", key, id, jsonPath); + } + const int k = std::stoi(key.substr(0, sep)); + const int l = std::stoi(key.substr(sep + 1)); + if (l < 1) { + LOGP(fatal, "Inextensional h key '{}' for sensor {}: l must be >= 1 (l = 0 is spanned by g) in {}", key, id, jsonPath); + } + sensor.inextensional.h[{k, l}] = val.get(); + } } - if (inex.contains("beta")) { - sensor.inextensional.beta = inex["beta"].get(); + // An "inextensional" block that yields no coefficients would silently + // produce a zero displacement field, i.e. a misalignment study that + // looks identical to the ideal one. Most likely cause: a file still in + // the old Fourier schema ("modes"/"alpha"/"beta"). + const auto& parsed = sensor.inextensional; + if (parsed.f.empty() && parsed.g.empty() && parsed.h.empty()) { + LOGP(fatal, + "Sensor {}: 'inextensional' block in {} contains none of the expected " + "keys 'f', 'g', 'h' - no deformation would be applied. Keys present: {}", + id, jsonPath, [&inex] { std::string s; for (const auto& [k, v] : inex.items()) { s += (s.empty() ? "" : ", ") + k; } return s; }()); } } } @@ -104,23 +129,6 @@ MisalignmentShift evaluateLegendreShift(const SensorMisalignment& sensor, const 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 x = {-0.9602898564975363, -0.7966664774136267, -0.5255324099163290, -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, 0.7966664774136267, 0.9602898564975363}; - constexpr std::array 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; @@ -137,31 +145,47 @@ MisalignmentShift evaluateInextensionalShift(const SensorMisalignment& sensor, c } const double r = constants::radii[frame.layerID]; - const double phi = std::atan2(r * std::sin(frame.alpha), r * std::cos(frame.alpha)); - const double z = frame.z; + const double gloX = frame.x * std::cos(frame.alpha); + const double gloY = frame.x * std::sin(frame.alpha); + const auto [u, v] = computeUV(gloX, gloY, frame.z, frame.sensorID, r); + const double cPhi = phiScale(r); + const double zOverR = frame.z / r; const auto& inex = sensor.inextensional; - double uz = 0., uphi = 0., ur = 0.; - for (const auto& [n, coeffs] : inex.modes) { - const double a_n = coeffs[0], b_n = coeffs[1], c_n = coeffs[2], d_n = coeffs[3]; - const double sn = std::sin(n * phi); - const double cn = std::cos(n * phi); - const int n2 = n * n; - - const double fn = (a_n * cn) + (b_n * sn); - const double fpn = (-n * a_n * sn) + (n * b_n * cn); - const double fppn = (-n2 * a_n * cn) - (n2 * b_n * sn); - const double gn = (c_n * cn) + (d_n * sn); - const double gpn = (-n * c_n * sn) + (n * d_n * cn); - - uz += fn; - uphi += -(z / r) * fpn + gn; - ur += (z / r) * fppn - gpn; + int maxK = 0; + for (const auto& [k, val] : inex.f) { + maxK = std::max(maxK, k); + } + for (const auto& [k, val] : inex.g) { + maxK = std::max(maxK, k); + } + int maxKh = 0, maxL = 0; + for (const auto& [kl, val] : inex.h) { + maxKh = std::max(maxKh, kl.first); + maxL = std::max(maxL, kl.second); } - uz += inex.alpha * phi; - uphi += -(z / r) * inex.alpha + inex.beta * phi; - ur += -inex.beta; + const auto pu = legendrePols(std::max(maxK, maxKh), u); + const auto pu1 = legendrePolsD1(maxK, u); + const auto pu2 = legendrePolsD2(maxK, u); + + // u_z = f, u_phi = -(z/r) f' + g, u_r = (z/r) f'' - g' + double uz = 0., uphi = 0., ur = 0.; + for (const auto& [k, fk] : inex.f) { + uz += fk * pu[k]; + uphi += -zOverR * cPhi * fk * pu1[k]; + ur += zOverR * cPhi * cPhi * fk * pu2[k]; + } + for (const auto& [k, gk] : inex.g) { + uphi += gk * pu[k]; + ur += -cPhi * gk * pu1[k]; + } + if (!inex.h.empty()) { + const auto pv = legendrePols(maxL, v); + for (const auto& [kl, hkl] : inex.h) { + ur += hkl * pu[kl.first] * pv[kl.second]; + } + } shift.dy = -uphi + (slopes.dydx * ur); shift.dz = -uz + (slopes.dzdx * ur); diff --git a/Detectors/Upgrades/ITS3/study/macros/PlotMisalignment.C b/Detectors/Upgrades/ITS3/study/macros/PlotMisalignment.C index 3f2170776e1fd..9c4e836bdf043 100644 --- a/Detectors/Upgrades/ITS3/study/macros/PlotMisalignment.C +++ b/Detectors/Upgrades/ITS3/study/macros/PlotMisalignment.C @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -74,6 +75,14 @@ void processTree(TFile* f, const char* treeName) const int nPtBins = 35; const double ptLimits[nPtBins] = {0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2., 2.2, 2.5, 3., 4., 5., 6., 8., 10., 15., 20.}; const int yDCABins{1000}; + const int nQ2PtBins = 2 * (nPtBins - 1); + double q2PtEdges[nQ2PtBins + 1]; + for (int i = 0; i < nPtBins; ++i) { + q2PtEdges[i] = -1.0 / ptLimits[i]; + } + for (int i = 1; i < nPtBins; ++i) { + q2PtEdges[nPtBins - 1 + i] = 1.0 / ptLimits[nPtBins - 1 - i]; + }; const float yDCARange{500}; // --- book histograms --- @@ -88,10 +97,14 @@ void processTree(TFile* f, const char* treeName) TH1F* hPull[kNPar][kNPtBins][kNLay]; // DCAxy TH2F* hDCAxyVsPt = new TH2F(Form("%s_hDCAxyVsPt", treeName), ";#it{p}_{T,MC} (GeV/#it{c});DCA_{#it{xy}} (#mum);entries", nPtBins - 1, ptLimits, yDCABins, -yDCARange, yDCARange); + TH2F* hDCAxyVsQ2Pt = new TH2F(Form("%s_hDCAxyVsQ2Pt", treeName), ";#it{q}/#it{p}_{T} (#it{c}/GeV);DCA_{#it{xy}} (#mum);entries", nQ2PtBins, q2PtEdges, yDCABins, -yDCARange, yDCARange); TH2F* hDCAxyVsPhi = new TH2F(Form("%s_hDCAxyVsPhi", treeName), ";#phi (rad);DCA_{#it{xy}} (#mum);entries", 100, 0, 2 * TMath::Pi(), yDCABins, -yDCARange, yDCARange); + TH2F* hDCAxyVsEta = new TH2F(Form("%s_hDCAxyVsEta", treeName), ";$eta;DCA_{#it{xy}} (#mum);entries", 200, -4, 4, yDCABins, -yDCARange, yDCARange); // DCAz TH2F* hDCAzVsPt = new TH2F(Form("%s_hDCAzVsPt", treeName), ";#it{p}_{T,MC} (GeV/#it{c});DCA_{#it{z}} (#mum);entries", nPtBins - 1, ptLimits, yDCABins, -yDCARange, yDCARange); TH2F* hDCAzVsPhi = new TH2F(Form("%s_hDCAzVsPhi", treeName), ";#phi (rad);DCA_{#it{z}} (#mum);entries", 100, 0, 2 * TMath::Pi(), yDCABins, -yDCARange, yDCARange); + TH2F* hDCAzVsEta = new TH2F(Form("%s_hDCAzVsEta", treeName), ";#eta;DCA_{#it{z}} (#mum);entries", 200, -4, 4, yDCABins, -yDCARange, yDCARange); + TH2F* hDCAzVsEtaAll = new TH2F(Form("%s_hDCAzVsEtaAll", treeName), ";#eta;DCA_{#it{z}} (#mum);entries", 200, -4, 4, yDCABins, -yDCARange, yDCARange); for (int ipt = 0; ipt < kNPtBins; ipt++) { for (int ilay = 0; ilay < kNLay; ilay++) { @@ -122,8 +135,8 @@ void processTree(TFile* f, const char* treeName) const Long64_t nEntries = tree->GetEntries(); for (Long64_t i = 0; i < nEntries; i++) { tree->GetEntry(i); - if (i % 100000 == 0) { - std::cout << "Progress: " << i << "/" << nEntries << " (" << (100.0 * i / nEntries) << "%)" << std::endl; + if (i % 1'000'000 == 0) { + std::cout << "Progress: " << i << "/" << nEntries << " (" << (100.0 * (double)i / (double)nEntries) << "%)" << '\n'; } int ilay = lay + 1; @@ -132,11 +145,19 @@ void processTree(TFile* f, const char* treeName) float dZum = dZ * 10000.f; if (lay == -1) { - hDCAxyVsPt->Fill(pt, dcaXY * 10000.); - hDCAzVsPt->Fill(pt, dcaZ * 10000.); + if (std::abs(eta) <= 1.0) { + hDCAxyVsPt->Fill(pt, dcaXY * 10000.); + hDCAxyVsQ2Pt->Fill(trk->getQ2Pt(), dcaXY * 10000.); + hDCAzVsPt->Fill(pt, dcaZ * 10000.); + hDCAzVsEtaAll->Fill(eta, dcaZ * 10000.); + if (pt >= 1.0 && pt <= 2.0) { + hDCAxyVsPhi->Fill(phi, dcaXY * 10000.); + hDCAzVsPhi->Fill(phi, dcaZ * 10000.); + } + } if (pt >= 1.0 && pt <= 2.0) { - hDCAxyVsPhi->Fill(phi, dcaXY * 10000.); - hDCAzVsPhi->Fill(phi, dcaZ * 10000.); + hDCAxyVsEta->Fill(eta, dcaXY * 10000.); + hDCAzVsEta->Fill(eta, dcaZ * 10000.); } } @@ -242,9 +263,13 @@ void processTree(TFile* f, const char* treeName) // write file out auto oFile = TFile::Open(Form("plotMisalignment_%s.root", treeName), "RECREATE"); hDCAxyVsPt->Write(); + hDCAxyVsQ2Pt->Write(); hDCAzVsPt->Write(); hDCAxyVsPhi->Write(); hDCAzVsPhi->Write(); + hDCAxyVsEta->Write(); + hDCAzVsEta->Write(); + hDCAzVsEtaAll->Write(); for (int ipt = 0; ipt < kNPtBins; ipt++) { for (int ilay = 0; ilay < kNLay; ilay++) { for (int iv = 0; iv < kNVar; iv++) { diff --git a/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx b/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx index c4e887ae315c3..04ba817fb5187 100644 --- a/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx +++ b/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx @@ -880,11 +880,17 @@ void TrackingStudySpec::doResidStudy() std::array cl; std::array clArr{nullptr}; dataformats::VertexBase pv; - float ip[2]; + float ip[2]{0.f, 0.f}; + // The MC vertex always defines the reference point of the impact + // parameters. Whether it *also* enters the track fit is a separate + // question, controlled by addPVAsCluster: constraining the fit to the same + // vertex the DCA is measured against compresses the DCA distribution, so + // for pointing-resolution studies the constraint must be switched off + // while the reference point stays defined. + const auto& eve = mMCReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); + pv.setXYZ(eve.GetX(), eve.GetY(), eve.GetZ()); if (mParams->addPVAsCluster) { - const auto& eve = mMCReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); auto trFitOut = iTrack.getParamIn(); - pv.setXYZ(eve.GetX(), eve.GetY(), eve.GetZ()); if (!prop->propagateToDCA(pv, trFitOut, bz, base::Propagator::MAX_STEP, mParams->CorrType)) { return; } @@ -927,9 +933,9 @@ void TrackingStudySpec::doResidStudy() } auto phi = i == 0 ? tInt.getPhi() : tInt.getPhiPos(); o2::math_utils::bringTo02Pi(phi); - if (clArr[0]) { - getImpactParams(tInt, pv, ip, bz); - } + // pv is always set, so the DCA no longer depends on the PV being + // part of the fit (it used to be gated on clArr[0]) + getImpactParams(tInt, pv, ip, bz); (*mDBGOut) << "res" << "dYInt=" << clArr[i]->getY() - tInt.getY() << "dZInt=" << clArr[i]->getZ() - tInt.getZ() @@ -979,54 +985,87 @@ void TrackingStudySpec::doMisalignmentStudy() int goodRefit{0}, notPassedSel{0}, fitFail{0}, fitFailMis{0}; o2::dataformats::VertexBase pv; - float ip[2]; + float ip[2]{0.f, 0.f}; float chi2{0}; auto writeTree = [&](const char* treeName, const std::array& clArr, const std::array& extrapOut, const std::array& extrapInw, const o2::MCCompLabel& lbl) { + const auto mcTrk = mMCReader.getTrack(lbl); + if (!mcTrk) { + return; + } + TParticlePDG* pPDG = TDatabasePDG::Instance()->GetParticle(mcTrk->GetPdgCode()); + if (!pPDG) { + return; + } + const std::array xyz{(float)mcTrk->GetStartVertexCoordinatesX(), (float)mcTrk->GetStartVertexCoordinatesY(), (float)mcTrk->GetStartVertexCoordinatesZ()}; + const std::array pxyz{(float)mcTrk->GetStartVertexMomentumX(), (float)mcTrk->GetStartVertexMomentumY(), (float)mcTrk->GetStartVertexMomentumZ()}; + const int mcCharge = TMath::Nint(pPDG->Charge() / 3); + + auto writeRow = [&](const o2::track::TrackParCov& tInt, float dY, float dZ, int layOut) { + // MC truth at same (alpha, x) + o2::track::TrackPar mcTrkAtX(xyz, pxyz, mcCharge, false); + if (!mcTrkAtX.rotate(tInt.getAlpha()) || !prop->PropagateToXBxByBz(mcTrkAtX, tInt.getX())) { + return; + } + auto phi = (layOut < 0) ? tInt.getPhi() : tInt.getPhiPos(); + o2::math_utils::bringTo02Pi(phi); + // pv is always set, so the DCA no longer depends on the PV being part of + // the fit (this used to be gated on clArr[0], which left dcaXY/dcaZ + // holding the previous track's values when the PV was disabled) + getImpactParams(tInt, pv, ip, prop->getNominalBz()); + (*mDBGOut) << treeName + << "trk=" << tInt + << "mcTrk=" << mcTrkAtX + << "chi2=" << chi2 + << "dY=" << dY + << "dZ=" << dZ + << "dcaXY=" << ip[0] + << "dcaZ=" << ip[1] + << "phi=" << phi + << "eta=" << tInt.getEta() + << "lay=" << layOut + << "\n"; + }; + + // Track-level row, lay = -1: the impact parameters with respect to the MC + // vertex. With addPVAsCluster the PV pseudo-cluster occupies slot 0 and + // supplies this row. Without it slot 0 is empty, so the row is emitted here + // from the innermost track state; otherwise every DCA histogram downstream + // (all of which select lay == -1) would come out empty. + // + // The interpolation at slot i deliberately excludes the cluster of slot i + // (that is what makes dY/dZ unbiased residuals), so it must be updated with + // that cluster before the impact parameters are taken. Without the update a + // misalignment confined to the innermost layer is invisible here by + // construction: the state would be built from the outer layers alone and + // would be bit-identical to the ideal one. + if (!clArr[0]) { + for (int i = 1; i <= 7; ++i) { + if (!clArr[i]) { + continue; + } + auto tFull = align::interpolateTrackParCov(extrapInw[i], extrapOut[i]); + if (!tFull.isValid() || !tFull.update(*clArr[i])) { + continue; + } + writeRow(tFull, 0.f, 0.f, -1); + break; + } + } + for (int i = 0; i <= 7; i++) { if (!clArr[i]) { continue; } // interpolated result - auto tInt = align::interpolateTrackParCov(extrapInw[i], extrapOut[i]); + const auto tInt = align::interpolateTrackParCov(extrapInw[i], extrapOut[i]); if (!tInt.isValid()) { continue; } - float dY = clArr[i]->getY() - tInt.getY(); - float dZ = clArr[i]->getZ() - tInt.getZ(); - // MC truth at same (alpha, x) - o2::track::TrackPar mcTrkAtX; - const auto mcTrk = mMCReader.getTrack(lbl); - if (mcTrk) { - std::array xyz{(float)mcTrk->GetStartVertexCoordinatesX(), (float)mcTrk->GetStartVertexCoordinatesY(), (float)mcTrk->GetStartVertexCoordinatesZ()}; - std::array pxyz{(float)mcTrk->GetStartVertexMomentumX(), (float)mcTrk->GetStartVertexMomentumY(), (float)mcTrk->GetStartVertexMomentumZ()}; - TParticlePDG* pPDG = TDatabasePDG::Instance()->GetParticle(mcTrk->GetPdgCode()); - if (pPDG) { - mcTrkAtX = o2::track::TrackPar(xyz, pxyz, TMath::Nint(pPDG->Charge() / 3), false); - if (mcTrkAtX.rotate(tInt.getAlpha()) && prop->PropagateToXBxByBz(mcTrkAtX, tInt.getX())) { - auto phi = i == 0 ? tInt.getPhi() : tInt.getPhiPos(); - o2::math_utils::bringTo02Pi(phi); - if (clArr[0]) { - getImpactParams(tInt, pv, ip, prop->getNominalBz()); - } - (*mDBGOut) << treeName - << "trk=" << tInt - << "mcTrk=" << mcTrkAtX - << "chi2=" << chi2 - << "dY=" << dY - << "dZ=" << dZ - << "dcaXY=" << ip[0] - << "dcaZ=" << ip[1] - << "phi=" << phi - << "eta=" << tInt.getEta() - << "lay=" << i - 1 - << "\n"; - } - } - } + writeRow(tInt, clArr[i]->getY() - tInt.getY(), clArr[i]->getZ() - tInt.getZ(), i - 1); } }; @@ -1052,12 +1091,14 @@ void TrackingStudySpec::doMisalignmentStudy() // ideal clusters std::array cl; std::array clArr{nullptr}; + // see the note in doRefits(): the MC vertex always defines the DCA + // reference point, addPVAsCluster only decides whether it constrains the fit + const auto& eve = mMCReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); + pv.setXYZ(eve.GetX(), eve.GetY(), eve.GetZ()); if (mParams->addPVAsCluster) { - const auto& eve = mMCReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); auto trFitOut = iTrack.getParamIn(); - pv.setXYZ(eve.GetX(), eve.GetY(), eve.GetZ()); if (!prop->propagateToDCA(pv, trFitOut, prop->getNominalBz(), base::Propagator::MAX_STEP, mParams->CorrType)) { - return; + continue; // skip this track, not the rest of the study } pv.setSigmaX(20e-4f); pv.setSigmaY(20e-4f);