diff --git a/Detectors/Upgrades/ITS3/alignment/CMakeLists.txt b/Detectors/Upgrades/ITS3/alignment/CMakeLists.txt index 0bc8080c7a1b8..e04dfcbb43963 100644 --- a/Detectors/Upgrades/ITS3/alignment/CMakeLists.txt +++ b/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 diff --git a/Detectors/Upgrades/ITS3/alignment/README.md b/Detectors/Upgrades/ITS3/alignment/README.md index 62633d1d7d313..80213eb4e03b1 100644 --- a/Detectors/Upgrades/ITS3/alignment/README.md +++ b/Detectors/Upgrades/ITS3/alignment/README.md @@ -27,4 +27,38 @@ dofSet.json: "calib": { "type": "legendre", "order": 1, "fix": [0, 2] } } ] -}``` +} +``` + + +## In-existensional modes +```json +{ + "defaults": { "rigidBody": "fixed" }, + "rules": [ + { + "match": "ITS3Layer1/ITS3CarbonForm0", + "calib": { + "type": "inextensional", + "order": 2, + "free": ["a_2", "b_2", "c_2", "d_2", "alpha", "beta"] + } + } + ] +} +``` + +```json +[ + { + "id": 2, + "inextensional": { + "modes": { + "2": [0.0008, -0.0005, 0.0006, -0.0007] + }, + "alpha": 0.0004, + "beta": -0.0003 + } + } +] +``` diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentDOF.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentDOF.h new file mode 100644 index 0000000000000..3fed9decbd6e7 --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentDOF.h @@ -0,0 +1,176 @@ +// 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_DOF_H +#define O2_ITS3_ALIGNMENT_DOF_H + +#include +#include +#include +#include +#include +#include + +#include + +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 +{ + public: + enum class Type : uint8_t { + RigidBody, + Legendre, + Inextensional + }; + virtual ~DOFSet() = default; + virtual Type type() const = 0; + int nDOFs() const { return static_cast(mFree.size()); } + virtual std::string dofName(int idx) const = 0; + virtual void fillDerivatives(const DerivativeContext& ctx, Eigen::Ref 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); } + int nFreeDOFs() const + { + int n = 0; + for (bool f : mFree) { + n += f; + } + return n; + } + + protected: + DOFSet(int n) : mFree(n, true) {} + std::vector mFree; +}; + +// Rigid body set +class RigidBodyDOFSet final : public DOFSet +{ + public: + // indices for rigid body parameters in LOC frame + enum RigidBodyDOF : uint8_t { + TX = 0, + TY, + TZ, + RX, + RY, + RZ, + NDOF, + }; + static constexpr const char* RigidBodyDOFNames[RigidBodyDOF::NDOF] = {"TX", "TY", "TZ", "RX", "RY", "RZ"}; + + RigidBodyDOFSet() : DOFSet(NDOF) {} + // mask: bitmask of free DOFs (bit i = DOF i is free) + explicit RigidBodyDOFSet(uint8_t mask) : DOFSet(NDOF) + { + for (int i = 0; i < NDOF; ++i) { + mFree[i] = (mask >> i) & 1; + } + } + Type type() const override { return Type::RigidBody; } + std::string dofName(int idx) const override { return RigidBodyDOFNames[idx]; } + void fillDerivatives(const DerivativeContext& ctx, Eigen::Ref out) const override; + uint8_t mask() const + { + uint8_t m = 0; + for (int i = 0; i < NDOF; ++i) { + m |= (uint8_t(mFree[i]) << i); + } + return m; + } +}; + +// Legendre DOFs +// Describing radial misplacement +class LegendreDOFSet final : public DOFSet +{ + public: + explicit LegendreDOFSet(int order) : DOFSet((order + 1) * (order + 2) / 2), mOrder(order) {} + Type type() const override { return Type::Legendre; } + int order() const { return mOrder; } + std::string dofName(int idx) const override + { + int i = 0; + while ((i + 1) * (i + 2) / 2 <= idx) { + ++i; + } + int j = idx - (i * (i + 1) / 2); + return std::format("L({},{})", i, j); + } + void fillDerivatives(const DerivativeContext& ctx, Eigen::Ref out) const override; + + private: + 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 +class InextensionalDOFSet final : public DOFSet +{ + public: + explicit InextensionalDOFSet(int maxOrder) : DOFSet((4 * (maxOrder - 1)) + 2), mMaxOrder(maxOrder) + { + if (maxOrder < 2) { + // the rest is eq. to rigid body + throw std::invalid_argument("InextensionalDOFSet requires maxOrder >= 2"); + } + } + Type type() const override { return Type::Inextensional; } + int maxOrder() const { return mMaxOrder; } + + // 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); } + + // indices of the non-periodic modes + int alphaIdx() const { return nPeriodic(); } + int betaIdx() const { return nPeriodic() + 1; } + + std::string dofName(int idx) const override + { + if (idx == alphaIdx()) { + return "alpha"; + } + if (idx == betaIdx()) { + return "beta"; + } + int n = (idx / 4) + 2; + int sub = idx % 4; + static constexpr const char* subNames[] = {"a", "b", "c", "d"}; + return std::format("{}_{}", subNames[sub], n); + } + void fillDerivatives(const DerivativeContext& ctx, Eigen::Ref out) const override; + + private: + int mMaxOrder; +}; + +#endif diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentHierarchy.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentHierarchy.h index 04b8157084d0a..ae8989deec21b 100644 --- a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentHierarchy.h +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentHierarchy.h @@ -13,14 +13,11 @@ #define O2_ITS3_ALIGNMENT_HIERARCHY_H #include -#include -#include -#include #include #include #include #include -#include +#include #include #include @@ -28,179 +25,14 @@ #include #include +#include "ITS3Align/AlignmentLabel.h" +#include "ITS3Align/AlignmentDOF.h" + namespace o2::its3::align { using Matrix36 = Eigen::Matrix; using Matrix66 = Eigen::Matrix; -// indices for rigid body parameters in LOC frame -enum RigidBodyDOF : uint8_t { - TX = 0, - TY, - TZ, - RX, - RY, - RZ, - NDOF, -}; -static constexpr const char* RigidBodyDOFNames[RigidBodyDOF::NDOF] = {"TX", "TY", "TZ", "RX", "RY", "RZ"}; - -// 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 DOFSet -{ - public: - enum class Type : uint8_t { RigidBody, - Legendre }; - virtual ~DOFSet() = default; - virtual Type type() const = 0; - int nDOFs() const { return static_cast(mFree.size()); } - virtual std::string dofName(int idx) 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); } - int nFreeDOFs() const - { - int n = 0; - for (bool f : mFree) { - n += f; - } - return n; - } - - protected: - DOFSet(int n) : mFree(n, true) {} - std::vector mFree; -}; - -class RigidBodyDOFSet final : public DOFSet -{ - public: - static constexpr int NDOF = RigidBodyDOF::NDOF; - RigidBodyDOFSet() : DOFSet(NDOF) {} - // mask: bitmask of free DOFs (bit i = DOF i is free) - explicit RigidBodyDOFSet(uint8_t mask) : DOFSet(NDOF) - { - for (int i = 0; i < NDOF; ++i) { - mFree[i] = (mask >> i) & 1; - } - } - Type type() const override { return Type::RigidBody; } - std::string dofName(int idx) const override { return RigidBodyDOFNames[idx]; } - uint8_t mask() const - { - uint8_t m = 0; - for (int i = 0; i < NDOF; ++i) { - m |= (uint8_t(mFree[i]) << i); - } - return m; - } -}; - -class LegendreDOFSet final : public DOFSet -{ - public: - explicit LegendreDOFSet(int order) : DOFSet((order + 1) * (order + 2) / 2), mOrder(order) {} - Type type() const override { return Type::Legendre; } - int order() const { return mOrder; } - std::string dofName(int idx) const override - { - int i = 0; - while ((i + 1) * (i + 2) / 2 <= idx) { - ++i; - } - int j = idx - (i * (i + 1) / 2); - return std::format("L({},{})", i, j); - } - - private: - int mOrder; -}; - -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) - 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 - static constexpr int ID_BITS = 22; // bits 6-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! - static constexpr T bitMask(int b) noexcept - { - return (T(1) << b) - T(1); - } - static constexpr int DOF_SHIFT = 0; - static constexpr T DOF_MAX = (T(1) << DOF_BITS) - T(1); - static constexpr T DOF_MASK = DOF_MAX << DOF_SHIFT; - static constexpr int CALIB_SHIFT = DOF_BITS; - static constexpr T CALIB_MAX = (T(1) << CALIB_BITS) - T(1); - static constexpr T CALIB_MASK = CALIB_MAX << CALIB_SHIFT; - static constexpr int ID_SHIFT = DOF_BITS + CALIB_BITS; - static constexpr T ID_MAX = (T(1) << ID_BITS) - T(1); - static constexpr T ID_MASK = ID_MAX << ID_SHIFT; - static constexpr int SENS_SHIFT = DOF_BITS + CALIB_BITS + ID_BITS; - static constexpr T SENS_MAX = (T(1) << SENS_BITS) - T(1); - static constexpr T SENS_MASK = SENS_MAX << SENS_SHIFT; - static constexpr int DET_SHIFT = DOF_BITS + CALIB_BITS + ID_BITS + SENS_BITS; - static constexpr T DET_MAX = (T(1) << DET_BITS) - T(1); - static constexpr T DET_MASK = DET_MAX << DET_SHIFT; - - 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)) - { - } - - /// produce the raw Millepede label for a given DOF index (rigid body: calib=0 in label) - constexpr T raw(T dof) const noexcept { return (mID & ~DOF_MASK) | ((dof & DOF_MAX) << DOF_SHIFT); } - constexpr int rawGBL(T dof) const noexcept { return static_cast(raw(dof)); } - - /// return a copy of this label with the CALIB bit set (for calibration DOFs on same volume) - GlobalLabel asCalib() const noexcept - { - GlobalLabel c{*this}; - c.mID |= (T(1) << CALIB_SHIFT); - return c; - } - - constexpr T id() const noexcept { return ((mID >> ID_SHIFT) & ID_MAX) - 1; } - constexpr T det() const noexcept { return (mID & DET_MASK) >> DET_SHIFT; } - constexpr bool sens() const noexcept { return (mID & SENS_MASK) >> SENS_SHIFT; } - constexpr bool calib() const noexcept { return (mID & CALIB_MASK) >> CALIB_SHIFT; } - - std::string asString() const - { - return std::format("Det:{} Id:{} Sens:{} Calib:{}", det(), id(), sens(), calib()); - } - - constexpr auto operator<=>(const GlobalLabel&) const noexcept = default; - - private: - T mID{0}; -}; - class HierarchyConstraint { public: @@ -220,8 +52,6 @@ class HierarchyConstraint std::vector mCoeff; // their coefficients }; -// --- AlignableVolume --- - class AlignableVolume { public: diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentLabel.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentLabel.h new file mode 100644 index 0000000000000..83495491b87e0 --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentLabel.h @@ -0,0 +1,87 @@ +// 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_LABEL_H +#define O2_ITS3_ALIGNMENT_LABEL_H + +#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) + 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 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! + static constexpr T bitMask(int b) noexcept + { + return (T(1) << b) - T(1); + } + static constexpr int DOF_SHIFT = 0; + static constexpr T DOF_MAX = (T(1) << DOF_BITS) - T(1); + static constexpr T DOF_MASK = DOF_MAX << DOF_SHIFT; + static constexpr int CALIB_SHIFT = DOF_BITS; + static constexpr T CALIB_MAX = (T(1) << CALIB_BITS) - T(1); + static constexpr T CALIB_MASK = CALIB_MAX << CALIB_SHIFT; + static constexpr int ID_SHIFT = DOF_BITS + CALIB_BITS; + static constexpr T ID_MAX = (T(1) << ID_BITS) - T(1); + static constexpr T ID_MASK = ID_MAX << ID_SHIFT; + static constexpr int SENS_SHIFT = DOF_BITS + CALIB_BITS + ID_BITS; + static constexpr T SENS_MAX = (T(1) << SENS_BITS) - T(1); + static constexpr T SENS_MASK = SENS_MAX << SENS_SHIFT; + static constexpr int DET_SHIFT = DOF_BITS + CALIB_BITS + ID_BITS + SENS_BITS; + static constexpr T DET_MAX = (T(1) << DET_BITS) - T(1); + static constexpr T DET_MASK = DET_MAX << DET_SHIFT; + + 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)) + { + } + + /// produce the raw Millepede label for a given DOF index (rigid body: calib=0 in label) + constexpr T raw(T dof) const noexcept { return (mID & ~DOF_MASK) | ((dof & DOF_MAX) << DOF_SHIFT); } + constexpr int rawGBL(T dof) const noexcept { return static_cast(raw(dof)); } + + /// return a copy of this label with the CALIB bit set (for calibration DOFs on same volume) + GlobalLabel asCalib() const noexcept + { + GlobalLabel c{*this}; + c.mID |= (T(1) << CALIB_SHIFT); + return c; + } + + constexpr T id() const noexcept { return ((mID >> ID_SHIFT) & ID_MAX) - 1; } + constexpr T det() const noexcept { return (mID & DET_MASK) >> DET_SHIFT; } + constexpr bool sens() const noexcept { return (mID & SENS_MASK) >> SENS_SHIFT; } + constexpr bool calib() const noexcept { return (mID & CALIB_MASK) >> CALIB_SHIFT; } + + std::string asString() const + { + return std::format("Det:{} Id:{} Sens:{} Calib:{}", det(), id(), sens(), calib()); + } + + constexpr auto operator<=>(const GlobalLabel&) const noexcept = default; + + private: + T mID{0}; +}; + +#endif diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentMath.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentMath.h new file mode 100644 index 0000000000000..9409648dca3e4 --- /dev/null +++ b/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 +#include + +namespace o2::its3::align +{ + +struct TrackSlopes { + double dydx{0.}; + double dzdx{0.}; +}; + +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); + +} // namespace o2::its3::align + +#endif diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentParams.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentParams.h index a7785a2c04e11..5a11066fd3c3b 100644 --- a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentParams.h +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentParams.h @@ -38,9 +38,10 @@ struct AlignmentParams : public o2::conf::ConfigurableParamHelper; diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/MisalignmentUtils.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/MisalignmentUtils.h new file mode 100644 index 0000000000000..457eccaeff4e6 --- /dev/null +++ b/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 +#include +#include +#include +#include + +#include "ITS3Align/AlignmentMath.h" +#include "MathUtils/LegendrePols.h" + +namespace o2::its3::align +{ + +struct InextensionalMisalignment { + std::map> 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 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 diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/TrackFit.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/TrackFit.h index 3f36705271c9b..4625776398c89 100644 --- a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/TrackFit.h +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/TrackFit.h @@ -76,10 +76,12 @@ o2::track::TrackParametrizationWithError interpolateTrackParCov( }; Mat55 cA = unpack(tA.getCov()); Mat55 cB = unpack(tB.getCov()); - Mat55 wA = cA.inverse(); - Mat55 wB = cB.inverse(); + Eigen::LLT lltA(cA), lltB(cB); + Mat55 wA = lltA.solve(Mat55::Identity()); + Mat55 wB = lltB.solve(Mat55::Identity()); Mat55 wTot = wA + wB; - Mat55 cTot = wTot.inverse(); + Eigen::LLT lltTot(wTot); + Mat55 cTot = lltTot.solve(Mat55::Identity()); Mat51 pA, pB; for (int i = 0; i < 5; ++i) { pA(i) = tA.getParam(i); diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentDOF.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentDOF.cxx new file mode 100644 index 0000000000000..d2a78dba791e6 --- /dev/null +++ b/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 +#include + +#include "ITS3Align/AlignmentMath.h" +#include "ITS3Base/SpecsV2.h" + +namespace +{ + +void validateDerivativeOutput(const DOFSet& dofSet, Eigen::Ref 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 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 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 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(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; +} diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx index 9170165a36a41..938c14c2c4759 100644 --- a/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx @@ -178,21 +178,23 @@ void AlignableVolume::writeParameters(std::ostream& os) const if (isRoot()) { os << "Parameter\n"; } - if (mRigidBody) { - for (int iDOF = 0; iDOF < mRigidBody->nDOFs(); ++iDOF) { - os << std::format("{:<10} {:>+15g} {:>+15g} ! {} {} ", - mLabel.raw(iDOF), 0.0, (mRigidBody->isFree(iDOF) ? 0.0 : -1.0), - (mRigidBody->isFree(iDOF) ? 'V' : 'F'), mRigidBody->dofName(iDOF)) - << mSymName << '\n'; + if (!mIsPseudo) { + if (mRigidBody) { + for (int iDOF = 0; iDOF < mRigidBody->nDOFs(); ++iDOF) { + os << std::format("{:<10} {:>+15g} {:>+15g} ! {} {} ", + mLabel.raw(iDOF), 0.0, (mRigidBody->isFree(iDOF) ? 0.0 : -1.0), + (mRigidBody->isFree(iDOF) ? 'V' : 'F'), mRigidBody->dofName(iDOF)) + << mSymName << '\n'; + } } - } - if (mCalib) { - auto calibLbl = mLabel.asCalib(); - for (int iDOF = 0; iDOF < mCalib->nDOFs(); ++iDOF) { - os << std::format("{:<10} {:>+15g} {:>+15g} ! {} {} ", - calibLbl.raw(iDOF), 0.0, (mCalib->isFree(iDOF) ? 0.0 : -1.0), - (mCalib->isFree(iDOF) ? 'V' : 'F'), mCalib->dofName(iDOF)) - << mSymName << '\n'; + if (mCalib) { + auto calibLbl = mLabel.asCalib(); + for (int iDOF = 0; iDOF < mCalib->nDOFs(); ++iDOF) { + os << std::format("{:<10} {:>+15g} {:>+15g} ! {} {:<5} ", + calibLbl.raw(iDOF), 0.0, (mCalib->isFree(iDOF) ? 0.0 : -1.0), + (mCalib->isFree(iDOF) ? 'V' : 'F'), mCalib->dofName(iDOF)) + << mSymName << '\n'; + } } } for (const auto& c : mChildren) { @@ -266,6 +268,9 @@ void applyDOFConfig(AlignableVolume* root, const std::string& jsonPath) } root->traverse([&](AlignableVolume* vol) { + if (vol->isPseudo()) { + return; + } const std::string& sym = vol->getSymName(); for (const auto& rule : rules) { const auto pattern = rule["match"].get(); @@ -357,6 +362,41 @@ 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); + bool fixed = cal.value("fixed", false); + if (fixed) { + dofSet->setAllFree(false); + } + if (cal.contains("free")) { + dofSet->setAllFree(false); + for (const auto& item : cal["free"]) { + if (item.is_number_integer()) { + dofSet->setFree(item.get(), true); + } else if (item.is_string()) { + for (int k = 0; k < dofSet->nDOFs(); ++k) { + if (dofSet->dofName(k) == item.get()) { + dofSet->setFree(k, true); + } + } + } + } + } + if (cal.contains("fix")) { + for (const auto& item : cal["fix"]) { + if (item.is_number_integer()) { + dofSet->setFree(item.get(), false); + } else if (item.is_string()) { + for (int k = 0; k < dofSet->nDOFs(); ++k) { + if (dofSet->dofName(k) == item.get()) { + dofSet->setFree(k, false); + } + } + } + } + } + vol->setCalib(std::move(dofSet)); } } } @@ -398,6 +438,12 @@ void writeMillepedeResults(AlignableVolume* root, const std::string& milleResPat // indexed by sensorID std::map> injRB; std::map>> injMatrix; + struct InjInex { + std::map> modes; + double alpha{0.}; + double beta{0.}; + }; + std::map injInex; if (!injectedJsonPath.empty()) { std::ifstream injFile(injectedJsonPath); if (injFile.is_open()) { @@ -410,6 +456,22 @@ void writeMillepedeResults(AlignableVolume* root, const std::string& milleResPat if (item.contains("matrix")) { injMatrix[id] = item["matrix"].get>>(); } + 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("alpha")) { + ii.alpha = inex["alpha"].get(); + } + if (inex.contains("beta")) { + ii.beta = inex["beta"].get(); + } + injInex[id] = ii; + } } LOGP(info, "Loaded injected misalignment for {} sensors", injData.size()); } else { @@ -468,6 +530,43 @@ void writeMillepedeResults(AlignableVolume* root, const std::string& milleResPat matrix.push_back(row); } entry["matrix"] = matrix; + } 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{}; + + 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; + } + 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; + + // beta + uint32_t rawBeta = calibLbl.raw(inexSet->betaIdx()); + auto itB = labelToValue.find(rawBeta); + inexEntry["beta"] = (itB != labelToValue.end() ? itB->second : 0.0) - inj.beta; + + entry["inextensional"] = inexEntry; } if (write) { output.push_back(entry); diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentMath.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentMath.cxx new file mode 100644 index 0000000000000..52e9c03540d4c --- /dev/null +++ b/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 + +#include + +#include "ITS3Base/SpecsV2.h" +#include "MathUtils/Utils.h" + +namespace o2::its3::align +{ + +std::pair 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 legendrePols(int order, double x) +{ + std::vector 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 diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentSpec.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentSpec.cxx index d381abc6aa567..72f968bdbf338 100644 --- a/Detectors/Upgrades/ITS3/alignment/src/AlignmentSpec.cxx +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentSpec.cxx @@ -10,8 +10,9 @@ // or submit itself to any jurisdiction. #include -#include #include +#include +#include #ifdef WITH_OPENMP #include @@ -43,12 +44,13 @@ #include "ITStracking/IOUtils.h" #include "ITS3Reconstruction/IOUtils.h" #include "ITS3Align/TrackFit.h" +#include "ITS3Align/AlignmentMath.h" #include "ITS3Align/AlignmentSpec.h" #include "ITS3Align/AlignmentParams.h" #include "ITS3Align/AlignmentTypes.h" #include "ITS3Align/AlignmentHierarchy.h" +#include "ITS3Align/MisalignmentUtils.h" #include "ITS3Align/AlignmentSensors.h" -#include "MathUtils/LegendrePols.h" namespace o2::its3::align { @@ -63,30 +65,29 @@ using TrackD = o2::track::TrackParCovD; namespace { -// compute normalized (u,v) in [-1,1] from global position on a sensor -std::pair computeUV(double gloX, double gloY, double gloZ, int sensorID, double radius) +DerivativeContext makeDerivativeContext(const FrameInfoExt& frame, const TrackD& trk) { - 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}; + const auto slopes = computeTrackSlopes(trk.getSnp(), trk.getTgl()); + const bool isITS3 = constants::detID::isDetITS3(frame.sens); + return {.sensorID = isITS3 ? constants::detID::getSensorID(frame.sens) : -1, + .layerID = isITS3 ? constants::detID::getDetID2Layer(frame.sens) : -1, + .measX = frame.x, + .measAlpha = frame.alpha, + .measZ = frame.positionTrackingFrame[1], + .trkY = trk.getY(), + .trkZ = trk.getZ(), + .snp = trk.getSnp(), + .tgl = trk.getTgl(), + .dydx = slopes.dydx, + .dzdx = slopes.dzdx}; } -// evaluate Legendre polynomials P_0(x) through P_order(x) via recurrence -std::vector legendrePols(int order, double x) +Matrix36 getRigidBodyBaseDerivatives(const DerivativeContext& ctx) { - std::vector 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; + static const RigidBodyDOFSet sRigidBodyBasis; + Eigen::MatrixXd dyn(3, sRigidBodyBasis.nDOFs()); + sRigidBodyBasis.fillDerivatives(ctx, dyn); + return dyn; } } // namespace @@ -152,8 +153,8 @@ class AlignmentSpec final : public Task GTrackID::mask_t mTracksSrc; int mNThreads{1}; const AlignmentParams* mParams{nullptr}; - std::array mDeformations; // one per sensorID (0-5) - std::array, 6> mRigidBodyParams; // (dx,dy,dz,rx,ry,rz) in LOC per sensorID + MisalignmentModel mMisalignment; + std::array, 6> mRigidBodyParams; // (dx,dy,dz,rx,ry,rz) in LOC per sensorID }; void AlignmentSpec::init(InitContext& ic) @@ -323,7 +324,8 @@ void AlignmentSpec::process() // this is the derivative in TRK but we want to align in LOC // so dr/da_(LOC) = dr/da_(TRK) * da_(TRK)/da_(LOC) const auto* tileVol = mChip2Hiearchy.at(lbl); - Matrix36 der = getRigidBodyDerivatives(wTrk); + const auto derCtx = makeDerivativeContext(frame, wTrk); + Matrix36 der = getRigidBodyBaseDerivatives(derCtx); // count rigid body columns: only volumes with real DOFs (not DOFPseudo) int nColRB{0}; @@ -375,39 +377,16 @@ void AlignmentSpec::process() } } - // 3) calibration derivatives (e.g. Legendre for ITS3 sensors, apply directly on the whole sensor, not on inidividual tiles) - if (calibSet && calibSet->type() == DOFSet::Type::Legendre) { - const auto* legSet = static_cast(calibSet); - const int N = legSet->order(); - const int sensorID = constants::detID::getSensorID(frame.sens); - const int layerID = constants::detID::getDetID2Layer(frame.sens); - - const double r = frame.x; - const double gX = r * std::cos(frame.alpha); - const double gY = r * std::sin(frame.alpha); - const double gZ = frame.positionTrackingFrame[1]; - auto [u, v] = computeUV(gX, gY, gZ, sensorID, constants::radii[layerID]); - - const double snp = wTrk.getSnp(); - const double tgl = wTrk.getTgl(); - const double csci = 1. / std::sqrt(1. - (snp * snp)); - const double dydx = snp * csci; - const double dzdx = tgl * csci; - - auto pu = legendrePols(N, u); - auto pv = legendrePols(N, v); - - int legIdx = 0; - const int legColStart = nColRB; - for (int i = 0; i <= N; ++i) { - for (int j = 0; j <= i; ++j) { - const double basis = pu[j] * pv[i - j]; - gLabels.push_back(sensorVol->getLabel().asCalib().rawGBL(legIdx)); - gDer(0, legColStart + legIdx) = dydx * basis; - gDer(1, legColStart + legIdx) = dzdx * basis; - ++legIdx; - } + // 3) calibration derivatives (apply directly on the whole sensor, not on individual tiles) + if (calibSet) { + const int nd = calibSet->nDOFs(); + Eigen::MatrixXd calDer(3, nd); + calibSet->fillDerivatives(derCtx, calDer); + for (int iDOF = 0; iDOF < nd; ++iDOF) { + gLabels.push_back(sensorVol->getLabel().asCalib().rawGBL(iDOF)); } + gDer.middleCols(curCol, nd) = calDer; + curCol += nd; } point.addGlobals(gLabels, gDer); } @@ -514,32 +493,24 @@ void AlignmentSpec::updateTimeDependentParams(ProcessingContext& pc) buildHierarchy(); - if (mParams->doMisalignmentLeg || mParams->doMisalignmentRB) { - TMatrixD null(1, 1); - null(0, 0) = 0; - for (int i = 0; i < 6; ++i) { - mDeformations[i] = o2::math_utils::Legendre2DPolynominal(null); - mRigidBodyParams[i].setZero(); + if (mParams->doMisalignmentLeg || mParams->doMisalignmentRB || mParams->doMisalignmentInex) { + mMisalignment = {}; + for (auto& rb : mRigidBodyParams) { + rb.setZero(); } if (!mParams->misAlgJson.empty()) { - using json = nlohmann::json; - std::ifstream f(mParams->misAlgJson); - auto data = json::parse(f); - for (const auto& item : data) { - int id = item["id"].get(); - if (mParams->doMisalignmentLeg && item.contains("matrix")) { - auto v = item["matrix"].get>>(); - TMatrixD m(v.size(), v[v.size() - 1].size()); - for (size_t r{0}; r < v.size(); ++r) { - for (size_t c{0}; c < v[r].size(); ++c) { - m(r, c) = v[r][c]; - } + mMisalignment = loadMisalignmentModel(mParams->misAlgJson); + if (mParams->doMisalignmentRB) { + using json = nlohmann::json; + std::ifstream f(mParams->misAlgJson); + auto data = json::parse(f); + for (const auto& item : data) { + int id = item["id"].get(); + if (!item.contains("rigidBody")) { + continue; } - mDeformations[id] = o2::math_utils::Legendre2DPolynominal(m); - } - if (mParams->doMisalignmentRB && item.contains("rigidBody")) { auto rb = item["rigidBody"].get>(); - for (int k = 0; k < 6 && k < (int)rb.size(); ++k) { + for (int k = 0; k < 6 && k < static_cast(rb.size()); ++k) { mRigidBodyParams[id](k) = rb[k]; } } @@ -848,6 +819,12 @@ bool AlignmentSpec::applyMisalignment(Eigen::Vector2d& res, const FrameInfoExt& const int sensorID = constants::detID::getSensorID(frame.sens); const int layerID = constants::detID::getDetID2Layer(frame.sens); + const MisalignmentFrame misFrame{ + .sensorID = sensorID, + .layerID = layerID, + .x = frame.x, + .alpha = frame.alpha, + .z = frame.positionTrackingFrame[1]}; // --- Legendre deformation (non-rigid-body) --- if (mParams->doMisalignmentLeg && mIsITS3 && mUseMC) { @@ -866,36 +843,18 @@ bool AlignmentSpec::applyMisalignment(Eigen::Vector2d& res, const FrameInfoExt& } o2::track::TrackParD mcPar(xyz, pxyz, TMath::Nint(pPDG->Charge() / 3), false); - const double r = frame.x; - const double gloX = r * std::cos(frame.alpha); - const double gloY = r * std::sin(frame.alpha); - const double gloZ = frame.positionTrackingFrame[1]; - auto [u, v] = computeUV(gloX, gloY, gloZ, sensorID, constants::radii[layerID]); - const double h = mDeformations[sensorID](u, v); - auto mcAtCl = mcPar; if (!mcAtCl.rotate(frame.alpha) || !prop->PropagateToXBxByBz(mcAtCl, frame.x)) { return false; } - const double snp = mcAtCl.getSnp(); - const double tgl = mcAtCl.getTgl(); - const double csci = 1. / std::sqrt(1. - (snp * snp)); - const double dydx = snp * csci; - const double dzdx = tgl * csci; - const double dy = dydx * h; - const double dz = dzdx * h; - - const double newGloY = (r * std::sin(frame.alpha)) + (dy * std::cos(frame.alpha)); - const double newGloX = (r * std::cos(frame.alpha)) - (dy * std::sin(frame.alpha)); - const double newGloZ = gloZ + dz; - auto [uNew, vNew] = computeUV(newGloX, newGloY, newGloZ, sensorID, constants::radii[layerID]); - if (std::abs(uNew) > 1. || std::abs(vNew) > 1.) { + const auto shift = evaluateLegendreShift(mMisalignment[sensorID], misFrame, computeTrackSlopes(mcAtCl.getSnp(), mcAtCl.getTgl())); + if (!shift.accepted) { return false; } - res[0] += dy; - res[1] += dz; + res[0] += shift.dy; + res[1] += shift.dz; } // --- Rigid body misalignment --- @@ -910,7 +869,7 @@ bool AlignmentSpec::applyMisalignment(Eigen::Vector2d& res, const FrameInfoExt& const auto* tileVol = mChip2Hiearchy.at(lbl); // derivative in TRK frame (3x6: rows = dy, dz, dsnp) - Matrix36 der = getRigidBodyDerivatives(wTrk); + Matrix36 der = getRigidBodyBaseDerivatives(makeDerivativeContext(frame, wTrk)); // TRK -> tile LOC const double posTrk[3] = {frame.x, 0., 0.}; @@ -929,6 +888,26 @@ bool AlignmentSpec::applyMisalignment(Eigen::Vector2d& res, const FrameInfoExt& res[1] += shift[1]; // dz } + // --- In-extensional deformation --- + // displacement field u(phi,z) = (u_phi, u_z, u_r) + // dy = -u_phi + y' * u_r, dz = -u_z + z' * u_r + if (mParams->doMisalignmentInex) { + const auto shift = evaluateInextensionalShift(mMisalignment[sensorID], misFrame, computeTrackSlopes(wTrk.getSnp(), wTrk.getTgl())); + res[0] += shift.dy; + res[1] += shift.dz; + } + + if (mOutOpt[OutputOpt::MisRes]) { + (*mDBGOut) << "mis" + << "dy=" << res[0] + << "dz=" << res[1] + << "sens=" << sensorID + << "lay=" << layerID + << "z=" << frame.positionTrackingFrame[1] + << "phi=" << frame.alpha + << "\n"; + } + return true; } diff --git a/Detectors/Upgrades/ITS3/alignment/src/MisalignmentUtils.cxx b/Detectors/Upgrades/ITS3/alignment/src/MisalignmentUtils.cxx new file mode 100644 index 0000000000000..ee5198ee98e0c --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/src/MisalignmentUtils.cxx @@ -0,0 +1,151 @@ +// 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/MisalignmentUtils.h" + +#include +#include +#include +#include + +#include +#include + +#include "Framework/Logger.h" +#include "ITS3Base/SpecsV2.h" + +namespace o2::its3::align +{ + +bool MisalignmentModel::empty() const noexcept +{ + return std::all_of(sensors.begin(), sensors.end(), [](const auto& sensor) { return sensor.empty(); }); +} + +MisalignmentModel loadMisalignmentModel(const std::string& jsonPath) +{ + MisalignmentModel model; + if (jsonPath.empty()) { + return model; + } + + std::ifstream f(jsonPath); + if (!f.is_open()) { + LOGP(fatal, "Cannot open misalignment JSON file: {}", jsonPath); + } + + using json = nlohmann::json; + const auto data = json::parse(f); + for (const auto& item : data) { + const int id = item["id"].get(); + if (id < 0 || id >= static_cast(MisalignmentModel::NSensors)) { + LOGP(fatal, "Misalignment sensor id {} out of range [0, {}) in {}", id, MisalignmentModel::NSensors, jsonPath); + } + + auto& sensor = model[id]; + if (item.contains("matrix")) { + auto v = item["matrix"].get>>(); + if (v.empty()) { + LOGP(fatal, "Legendre matrix for sensor {} is empty in {}", id, jsonPath); + } + TMatrixD m(v.size(), v.back().size()); + for (std::size_t r{0}; r < v.size(); ++r) { + for (std::size_t c{0}; c < v[r].size(); ++c) { + m(r, c) = v[r][c]; + } + } + sensor.legendre = o2::math_utils::Legendre2DPolynominal(m); + sensor.hasLegendre = true; + } + 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>(); + } + } + if (inex.contains("alpha")) { + sensor.inextensional.alpha = inex["alpha"].get(); + } + if (inex.contains("beta")) { + sensor.inextensional.beta = inex["beta"].get(); + } + } + } + + return model; +} + +MisalignmentShift evaluateLegendreShift(const SensorMisalignment& sensor, const MisalignmentFrame& frame, const TrackSlopes& slopes) +{ + MisalignmentShift shift; + if (!sensor.hasLegendre) { + return shift; + } + + const double gloX = frame.x * std::cos(frame.alpha); + const double gloY = frame.x * std::sin(frame.alpha); + const double gloZ = frame.z; + auto [u, v] = computeUV(gloX, gloY, gloZ, frame.sensorID, constants::radii[frame.layerID]); + const double h = sensor.legendre(u, v); + + shift.dy = slopes.dydx * h; + shift.dz = slopes.dzdx * h; + + 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; + auto [uNew, vNew] = computeUV(newGloX, newGloY, newGloZ, frame.sensorID, constants::radii[frame.layerID]); + shift.accepted = std::abs(uNew) <= 1. && std::abs(vNew) <= 1.; + return shift; +} + +MisalignmentShift evaluateInextensionalShift(const SensorMisalignment& sensor, const MisalignmentFrame& frame, const TrackSlopes& slopes) +{ + MisalignmentShift shift; + if (!sensor.hasInextensional) { + return shift; + } + + 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 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; + } + + uz += inex.alpha * phi; + uphi += -(z / r) * inex.alpha + inex.beta * phi; + ur += -inex.beta; + + shift.dy = -uphi + (slopes.dydx * ur); + shift.dz = -uz + (slopes.dzdx * ur); + return shift; +} + +} // namespace o2::its3::align diff --git a/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx b/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx index 4ce2c79cb23f1..63a08d1c9c6ab 100644 --- a/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx +++ b/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx @@ -11,12 +11,10 @@ #include #include -#include #include #include #include -#include #include "CommonUtils/TreeStreamRedirector.h" #include "DataFormatsGlobalTracking/RecoContainer.h" @@ -38,6 +36,7 @@ #include "ITS3Reconstruction/IOUtils.h" #include "ITS3TrackingStudy/ITS3TrackingStudyParam.h" #include "ITS3TrackingStudy/ParticleInfoExt.h" +#include "ITS3Align/MisalignmentUtils.h" #include "ITS3Align/TrackFit.h" #include "ReconstructionDataFormats/DCA.h" #include "ReconstructionDataFormats/GlobalTrackID.h" @@ -47,7 +46,6 @@ #include "SimulationDataFormat/MCEventLabel.h" #include "SimulationDataFormat/MCUtils.h" #include "Steer/MCKinematicsReader.h" -#include "MathUtils/LegendrePols.h" #include "Framework/Logger.h" namespace o2::its3::study @@ -153,7 +151,7 @@ class TrackingStudySpec : public Task o2::steer::MCKinematicsReader mMCReader; // reader of MC information const o2::its3::TopologyDictionary* mITSDict = nullptr; // cluster patterns dictionary o2::globaltracking::RecoContainer mRecoData; - std::array mDeformations; // one per sensorID (0-5) + align::MisalignmentModel mMisalignment; }; void TrackingStudySpec::init(InitContext& ic) @@ -186,26 +184,9 @@ void TrackingStudySpec::updateTimeDependentParams(ProcessingContext& pc) o2::its::GeometryTGeo::Instance()->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G, o2::math_utils::TransformType::T2G)); mParams = &ITS3TrackingStudyParam::Instance(); if (mParams->doMisalignment) { - TMatrixD null(1, 1); - null(0, 0) = 0; - for (int i = 0; i < 6; ++i) { - mDeformations[i] = o2::math_utils::Legendre2DPolynominal(null); - } + mMisalignment = {}; if (!mParams->misAlgJson.empty()) { - using json = nlohmann::json; - std::ifstream f(mParams->misAlgJson); - auto data = json::parse(f); - for (const auto& item : data) { - int id = item["id"].get(); - auto v = item["matrix"].get>>(); - TMatrixD m(v.size(), v[v.size() - 1].size()); - for (size_t r{0}; r < v.size(); ++r) { - for (size_t c{0}; c < v[r].size(); ++c) { - m(r, c) = v[r][c]; - } - } - mDeformations[id] = o2::math_utils::Legendre2DPolynominal(m); - } + mMisalignment = align::loadMisalignmentModel(mParams->misAlgJson); } } } @@ -984,17 +965,6 @@ void TrackingStudySpec::doMisalignmentStudy() int goodRefit{0}, notPassedSel{0}, fitFail{0}, fitFailMis{0}; - // compute normalized (u,v) in [-1,1] from global position - auto computeUV = [](float gloX, float gloY, float gloZ, int sensorID, float radius) -> std::pair { - const bool isTop = sensorID % 2 == 0; - const double phi = o2::math_utils::to02Pi(std::atan2(gloY, gloX)); - const double phiBorder1 = o2::math_utils::to02Pi(((isTop ? 0. : 1.) * TMath::Pi()) + std::asin(constants::equatorialGap / 2. / radius)); - const double phiBorder2 = o2::math_utils::to02Pi(((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}; - }; - float chi2{0}; auto writeTree = [&](const char* treeName, const std::array& clArr, @@ -1102,8 +1072,8 @@ void TrackingStudySpec::doMisalignmentStudy() } writeTree("idealRes", clArr, extrapOut, extrapInw, lbl); - // Propagate MC truth to each cluster's (alpha, x) to get true track direction, - // then compute dy = dydx*h(u,v), dz = dzdx*h(u,v) - first Newton step. + // Propagate MC truth to each cluster's (alpha, x) to get true track direction. + // The shared misalignment evaluators then provide the tracking-frame dy/dz shift. const auto mcTrk = mMCReader.getTrack(lbl); if (!mcTrk) { continue; @@ -1132,14 +1102,7 @@ void TrackingStudySpec::doMisalignmentStudy() } const int sensorID = constants::detID::getSensorID(sens); const int layerID = constants::detID::getDetID2Layer(sens); - - // compute h(u,v) at this cluster - const float r = orig.getX(); - const float gloX = r * std::cos(orig.alpha); - const float gloY = r * std::sin(orig.alpha); - const float gloZ = orig.getZ(); - auto [u, v] = computeUV(gloX, gloY, gloZ, sensorID, constants::radii[layerID]); - const double h = mDeformations[sensorID](u, v); + const auto& sensorMis = mMisalignment[sensorID]; // propagate MC track to cluster's tracking frame to get true slopes auto mcAtCl = mcPar; @@ -1147,28 +1110,31 @@ void TrackingStudySpec::doMisalignmentStudy() clArrMis[1 + iLay] = nullptr; // can't compute slopes -> drop cluster continue; } - const float snp = mcAtCl.getSnp(); - const float tgl = mcAtCl.getTgl(); - const float csci = 1.f / std::sqrt(1.f - (snp * snp)); - const float dydx = snp * csci; - const float dzdx = tgl * csci; - const float dy = dydx * static_cast(h); - const float dz = dzdx * static_cast(h); - - // check if shifted position is still within sensor acceptance - const float newGloY = (r * std::sin(orig.alpha)) + (dy * std::cos(orig.alpha)); - const float newGloX = (r * std::cos(orig.alpha)) - (dy * std::sin(orig.alpha)); - const float newGloZ = gloZ + dz; - auto [uNew, vNew] = computeUV(newGloX, newGloY, newGloZ, sensorID, constants::radii[layerID]); - if (std::abs(uNew) > 1. || std::abs(vNew) > 1.) { - clArrMis[1 + iLay] = nullptr; // shifted outside acceptance - continue; + const align::MisalignmentFrame misFrame{ + .sensorID = sensorID, + .layerID = layerID, + .x = orig.getX(), + .alpha = orig.alpha, + .z = orig.getZ()}; + const auto slopes = align::computeTrackSlopes(mcAtCl.getSnp(), mcAtCl.getTgl()); + + align::MisalignmentShift totalShift; + if (sensorMis.hasLegendre) { + const auto shift = align::evaluateLegendreShift(sensorMis, misFrame, slopes); + if (!shift.accepted) { + clArrMis[1 + iLay] = nullptr; // shifted outside acceptance + continue; + } + totalShift += shift; + } + if (sensorMis.hasInextensional) { + totalShift += align::evaluateInextensionalShift(sensorMis, misFrame, slopes); } // create shifted copy: keep x=r (nominal), shift y and z misClArr[iLay] = orig; - misClArr[iLay].setY(orig.getY() + dy); - misClArr[iLay].setZ(orig.getZ() + dz); + misClArr[iLay].setY(orig.getY() + totalShift.dy); + misClArr[iLay].setZ(orig.getZ() + totalShift.dz); misClArr[iLay].setSigmaY2(orig.getSigmaY2() + (mParams->misAlgExtCY[sensorID] * mParams->misAlgExtCY[sensorID])); misClArr[iLay].setSigmaZ2(orig.getSigmaZ2() + (mParams->misAlgExtCZ[sensorID] * mParams->misAlgExtCZ[sensorID])); clArrMis[1 + iLay] = &misClArr[iLay];