ITS3: alignment

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

Author: f3sch

Detectors/Upgrades/ITS3/alignment/CMakeLists.txt

@@ -9,18 +9,42 @@
 # granted to it by virtue of its status as an Intergovernmental Organization
 # or submit itself to any jurisdiction.
 
+# add_compile_options(-O0 -g -fPIC -fno-omit-frame-pointer)
 o2_add_library(ITS3Align
+               TARGETVARNAME targetName
                SOURCES src/MisalignmentParameters.cxx
                        src/MisalignmentHits.cxx
                        src/MisalignmentManager.cxx
                        src/Deformations.cxx
+                       src/AlignmentHierarchy.cxx
+                       src/AlignmentSensors.cxx
+                       src/AlignmentParams.cxx
+                       src/AlignmentTypes.cxx
+                       src/AlignmentSpec.cxx
                PUBLIC_LINK_LIBRARIES O2::MathUtils
                                      O2::Steer
                                      O2::ITSBase
-                                     O2::ITSMFTSimulation)
+                                     O2::ITSMFTSimulation
+                                     O2::Framework
+                                     O2::GlobalTrackingWorkflowReaders
+                                     O2::GlobalTrackingWorkflowHelpers
+                                     O2::DataFormatsGlobalTracking
+                                     O2::DetectorsVertexing
+                                     GBL::GBL)
+if (OpenMP_CXX_FOUND)
+    target_compile_definitions(${targetName} PRIVATE WITH_OPENMP)
+    target_link_libraries(${targetName} PRIVATE OpenMP::OpenMP_CXX)
+endif()
 
 o2_target_root_dictionary(ITS3Align
                           HEADERS include/ITS3Align/MisalignmentParameters.h
                                   include/ITS3Align/MisalignmentHits.h
-                                  include/ITS3Align/MisalignmentHits.h
-                                  include/ITS3Align/Deformations.h)
+                                  include/ITS3Align/Deformations.h
+                                  include/ITS3Align/AlignmentParams.h
+                                  include/ITS3Align/AlignmentTypes.h)
+
+
+o2_add_executable(alignment-workflow
+                  SOURCES src/alignment-workflow.cxx
+                  COMPONENT_NAME its3
+                  PUBLIC_LINK_LIBRARIES O2::ITS3Align)

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

@@ -0,0 +1,258 @@
+// 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_HIERARCHY_H
+#define O2_ITS3_ALIGNMENT_HIERARCHY_H
+
+#include <memory>
+#include <compare>
+#include <type_traits>
+#include <map>
+#include <utility>
+#include <vector>
+#include <ostream>
+
+#include <Eigen/Dense>
+
+#include <TGeoMatrix.h>
+#include <TGeoPhysicalNode.h>
+
+namespace o2::its3::align
+{
+using Matrix36 = Eigen::Matrix<double, 3, 6>;
+using Matrix66 = Eigen::Matrix<double, 6, 6>;
+
+using RigidBodyDOFMask = uint8_t;
+template <typename... Args>
+constexpr RigidBodyDOFMask DOF_BIT(Args... bits)
+{
+  return (RigidBodyDOFMask{0} | ... | (RigidBodyDOFMask{1} << bits));
+}
+// DOFs are defined in LOC
+enum RigidBodyDOF : RigidBodyDOFMask {
+  TX = 0, // Translation along local X
+  TY,     // Translation along local Y
+  TZ,     // Translation along local Z
+  RX,     // Rotation around local X
+  RY,     // Rotation around local Y
+  RZ,     // Rotation around local Z
+  NDOF,
+};
+constexpr RigidBodyDOFMask RigidBodyDOFTra = DOF_BIT(TX, TY, TZ);
+constexpr RigidBodyDOFMask RigidBodyDOFRot = DOF_BIT(RX, RY, RZ);
+constexpr RigidBodyDOFMask RigidBodyDOFAll = RigidBodyDOFTra | RigidBodyDOFRot;
+constexpr RigidBodyDOFMask RigidBodyDOFNone = 0;
+constexpr RigidBodyDOFMask RigidBodyDOFPseudo = std::numeric_limits<RigidBodyDOFMask>::max(); // special value representing that the volume itself does not have any dofs but should not curtail the parent's ones
+static constexpr const char* RigidBodyDOFNames[RigidBodyDOF::NDOF] = {"TX", "TY", "TZ", "RX", "RY", "RZ"};
+inline bool hasRigidBodyDOF(RigidBodyDOFMask m, RigidBodyDOFMask d)
+{
+  return (m == RigidBodyDOFPseudo) || (m & DOF_BIT(d));
+}
+inline void enableRigidBodyDOF(RigidBodyDOFMask& m, RigidBodyDOFMask d)
+{
+  m |= DOF_BIT(d);
+}
+inline void disableRigidBodyDOF(RigidBodyDOFMask& m, RigidBodyDOFMask d)
+{
+  m &= ~DOF_BIT(d);
+}
+// 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 GlobalLabel
+{
+  // Millepede label is any positive integer [1....)
+ public:
+  using T = uint32_t;
+  static constexpr int DOF_BITS = 5;  // 0...4
+  static constexpr int ID_BITS = 23;  // 5...27
+  static constexpr int SENS_BITS = 1; // 28
+  static constexpr int TOTAL_BITS = sizeof(T) * 8;
+  static constexpr int DET_BITS = TOTAL_BITS - (DOF_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 ID_SHIFT = DOF_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 = ID_BITS + DOF_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 + 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) : mID(((((id + 1) & ID_MAX) << ID_SHIFT) | ((det & DET_MAX) << DET_SHIFT) | ((sens & SENS_MAX) << SENS_SHIFT)))
+  {
+  }
+
+  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<int>(raw(dof)); }
+  constexpr T id() const noexcept
+  {
+    return ((mID - 1) & ID_MASK) >> ID_SHIFT;
+  }
+  constexpr T det() const noexcept
+  {
+    return (mID & DET_MASK) >> DET_SHIFT;
+  }
+  constexpr bool sens() const noexcept
+  {
+    return (mID & SENS_MASK) >> SENS_SHIFT;
+  }
+
+  std::string asString() const
+  {
+    return std::format("Det:{} Id:{} Sens:{}", det(), id(), sens());
+  }
+
+  constexpr auto operator<=>(const GlobalLabel&) const noexcept = default;
+
+ private:
+  T mID{0};
+};
+
+// Rigid body constraints for the parents
+class HierarchyConstraint
+{
+ public:
+  HierarchyConstraint(std::string name, double value) : mName(std::move(name)), mValue(value) {}
+  void add(uint32_t lab, double coeff)
+  {
+    mLabels.push_back(lab);
+    mCoeff.push_back(coeff);
+  }
+  void write(std::ostream& os) const;
+  auto getSize() const noexcept { return mLabels.size(); }
+
+ private:
+  std::string mName;             // name of the constraint
+  double mValue{0.0};            // constraint value
+  std::vector<uint32_t> mLabels; // parameter labels
+  std::vector<double> mCoeff;    // their coefficients
+};
+
+class AlignableVolume
+{
+ public:
+  using Ptr = std::unique_ptr<AlignableVolume>;
+  using SensorMapping = std::map<GlobalLabel, AlignableVolume*>;
+
+  AlignableVolume(const AlignableVolume&) = default;
+  AlignableVolume(AlignableVolume&&) = delete;
+  AlignableVolume& operator=(const AlignableVolume&) = default;
+  AlignableVolume& operator=(AlignableVolume&&) = delete;
+  AlignableVolume(const char* symName, uint32_t label, uint32_t det, bool sens, RigidBodyDOFMask dof);
+  AlignableVolume(const char* symName, GlobalLabel label, RigidBodyDOFMask dof);
+  virtual ~AlignableVolume() = default;
+
+  void finalise(uint8_t level = 0);
+
+  // steering file output
+  void writeRigidBodyConstraints(std::ostream& os) const;
+  void writeParameters(std::ostream& os) const;
+  void writeTree(std::ostream& os, int indent = 0) const;
+
+  // tree-like
+  auto getLevel() const noexcept { return mLevel; }
+  bool isRoot() const noexcept { return mParent == nullptr; }
+  bool isLeaf() const noexcept { return mChildren.empty(); }
+  template <class T = AlignableVolume>
+    requires std::derived_from<T, AlignableVolume>
+  AlignableVolume* addChild(const char* symName, uint32_t label, uint32_t det, bool sens, RigidBodyDOFMask dof)
+  {
+    auto c = std::make_unique<T>(symName, label, det, sens, dof);
+    return setParent(std::move(c));
+  }
+  template <class T = AlignableVolume>
+    requires std::derived_from<T, AlignableVolume>
+  AlignableVolume* addChild(const char* symName, GlobalLabel lbl, RigidBodyDOFMask dof)
+  {
+    auto c = std::make_unique<T>(symName, lbl, dof);
+    return setParent(std::move(c));
+  }
+
+  // bfs traversal
+  void traverse(const std::function<void(AlignableVolume*)>& visitor)
+  {
+    visitor(this);
+    for (auto& c : mChildren) {
+      c->traverse(visitor);
+    }
+  }
+
+  std::string getSymName() const noexcept { return mSymName; }
+  GlobalLabel getLabel() const noexcept { return mLabel; }
+  AlignableVolume* getParent() const { return mParent; }
+  size_t getNChildren() const noexcept { return mChildren.size(); }
+  void setRigidBodyDOF(RigidBodyDOFMask m) noexcept { mDOF = m; }
+  RigidBodyDOFMask getRigidBodyDOF() const noexcept { return mDOF; }
+
+  // transformation matrices
+  virtual void defineMatrixL2G() {}
+  virtual void defineMatrixT2L() {}
+  virtual void computeJacobianL2T(const double* pos, Matrix66& jac) const {};
+  const TGeoHMatrix& getL2P() const { return mL2P; }
+  const TGeoHMatrix& getT2L() const { return mT2L; }
+  const Matrix66& getJL2P() const { return mJL2P; }
+  const Matrix66& getJP2L() const { return mJP2L; }
+
+ protected:
+  /// matrices
+  AlignableVolume* mParent{nullptr}; // parent
+  TGeoPNEntry* mPNE{nullptr};        // physical entry
+  TGeoPhysicalNode* mPN{nullptr};    // physical node
+  TGeoHMatrix mL2G;                  // (LOC) -> (GLO)
+  TGeoHMatrix mL2P;                  // (LOC) -> (PAR)
+  Matrix66 mJL2P;                    // jac (LOC) -> (PAR)
+  Matrix66 mJP2L;                    // jac (PAR) -> (LOC)
+  TGeoHMatrix mT2L;                  // (TRK) -> (LOC)
+
+ private:
+  double mPreSigma{0.0};                  // asigned pre-sigma
+  std::string mSymName;                   // unique symbolic name
+  GlobalLabel mLabel;                     // internal global idetenifier
+  uint8_t mLevel{0};                      // depth-in tree
+  RigidBodyDOFMask mDOF{RigidBodyDOFAll}; // allowed dofs
+
+  AlignableVolume* setParent(Ptr c)
+  {
+    c->mParent = this;
+    mChildren.push_back(std::move(c));
+    return mChildren.back().get();
+  }
+  std::vector<Ptr> mChildren; // children
+
+  void init();
+};
+
+} // namespace o2::its3::align
+
+#endif

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

@@ -0,0 +1,54 @@
+// Copyright 2019-2020 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 ALICEO2_ITS3_ALIGNMENTPARAMS_H_
+#define ALICEO2_ITS3_ALIGNMENTPARAMS_H_
+
+#include "CommonUtils/ConfigurableParam.h"
+#include "CommonUtils/ConfigurableParamHelper.h"
+#include "DetectorsBase/Propagator.h"
+
+namespace o2::its3::align
+{
+struct AlignmentParams : public o2::conf::ConfigurableParamHelper<AlignmentParams> {
+  // Track selection
+  float minPt = 1.f;         // minimum pt required
+  int minITSCls = 7;         // minimum number of ITS clusters
+  float maxITSChi2Ndf = 1.2; // maximum ITS track chi2
+
+  // propagation opt
+  double maxSnp = 0.85;
+  double maxStep = 2.0;
+  // o2::base::PropagatorD::MatCorrType matCorrType = o2::base::PropagatorD::MatCorrType::USEMatCorrTGeo;
+  o2::base::PropagatorD::MatCorrType matCorrType = o2::base::PropagatorD::MatCorrType::USEMatCorrLUT;
+
+  bool useStableRefit = true; // use input tracks as linearization point
+  float minMS = 1e-6f;        //  minimum scattering to account for
+  float maxChi2Ndf = 10;      // maximum Chi2/Ndf allowed for GBL fit
+
+  // Ridder options
+  int ridderMaxExtrap = 10;
+  double ridderRelIniStep[5] = {0.01, 0.01, 0.02, 0.02, 0.02};
+  double ridderMaxIniStep[5] = {0.1, 0.1, 0.05, 0.05, 0.05};
+  double ridderShrinkFac = 2.0;
+  double ridderEps = 1e-16;
+  double ridderMaxJacDiagTol = 0.1; // max tolerance of diagonal elements away from 1
+
+  // MillePede output
+  std::string milleBinFile = "mp2data.bin";
+  std::string milleConFile = "mp2con.txt";
+  std::string milleParamFile = "mp2param.txt";
+  std::string milleTreeFile = "mp2tree.txt";
+
+  O2ParamDef(AlignmentParams, "ITS3AlignmentParams");
+};
+} // namespace o2::its3::align
+
+#endif