something with ITS3

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

Author: f3sch

Detectors/ITSMFT/ITS/simulation/include/ITSSimulation/Detector.h

@@ -216,7 +216,7 @@ class Detector : public o2::base::DetImpl<Detector>
 
   /// Return Chip Volume UID
   /// \param id volume id
-  Int_t chipVolUID(Int_t id) const { return o2::base::GeometryManager::getSensID(o2::detectors::DetID::ITS, id); }
+  Int_t chipVolUID(Int_t id) const { return o2::base::GeometryManager::getSensID(o2::detectors::DetID(GetName()), id); }
 
   void EndOfEvent() override;
 

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

@@ -14,6 +14,7 @@
 
 #include "ITS3Align/MisalignmentParameters.h"
 #include "MathUtils/LegendrePols.h"
+#include "GPUCommonMath.h"
 
 #include <filesystem>
 
@@ -26,57 +27,31 @@ class Deformations
   // init deformations from the parameter file
   void init(const std::filesystem::path&);
 
-  double getDeformationX(unsigned int id, double u, double v) const { return getDeformation<0>(id, u, v); }
-  double getDeformationY(unsigned int id, double u, double v) const { return getDeformation<1>(id, u, v); }
-  double getDeformationZ(unsigned int id, double u, double v) const { return getDeformation<2>(id, u, v); }
-  double getDeformation(unsigned int id, unsigned int axis, double u, double v) const
+  std::tuple<double, double, double> getDeformation(unsigned int id, double radius, double phi, double z, double u, double v, double fac = 1.0) const
   {
-    if (axis == 0) {
-      return mLegX[id](u, v);
-    } else if (axis == 1) {
-      return mLegY[id](u, v);
-    } else {
-      return mLegZ[id](u, v);
-    }
+    // Calculate f_def(phi,z) = ((r+dr)*cos(phi), (r+dr)*sin(phi), z)^T + (dx, dy, dz)^T
+    const double dr = mLegendre[id](u, v);
+    const double drr = radius + dr * fac;
+    double sn, cs;
+    o2::gpu::GPUCommonMath::SinCosd(phi, sn, cs);
+    const auto& global = mParams.getGlobal(id);
+    return {drr * cs + global.getX() * fac, drr * sn + global.getY() * fac, z + global.getZ() * fac};
   }
-  std::array<double, 3> getDeformation(unsigned int id, double u, double v) const
-  {
-    return {getDeformation<0>(id, u, v),
-            getDeformation<1>(id, u, v),
-            getDeformation<2>(id, u, v)};
-  }
-  std::array<unsigned int, 3> getOrders(unsigned int id) const
-  {
-    return {mLegX[id].NOrder(), mLegY[id].NOrder(), mLegZ[id].NOrder()};
-  }
-  const o2::math_utils::Legendre2DPolynominal& getLegendre(unsigned int id, unsigned int axis) const
+
+  double getDeformation(unsigned int id, double u, double v)
   {
-    if (axis == 0) {
-      return mLegX[id];
-    } else if (axis == 1) {
-      return mLegY[id];
-    } else {
-      return mLegZ[id];
-    }
+    return mLegendre[id](u, v);
   }
 
+  const o2::math_utils::Legendre2DPolynominal& getLegendre(unsigned int id) { return mLegendre[id]; }
+
+  unsigned int getOrder(unsigned int id) { return mLegendre[id].NOrder(); }
+
  private:
-  template <int axis>
-  double getDeformation(unsigned int id, double u, double v) const
-  {
-    if constexpr (axis == 0) {
-      return mLegX[id](u, v);
-    } else if constexpr (axis == 1) {
-      return mLegY[id](u, v);
-    } else {
-      return mLegZ[id](u, v);
-    }
-  }
+  MisalignmentParameters mParams;
 
-  // 3 Legendre polynominals to model deformations in x,y,z; parameterized by normalized phi (u) and z (v) coordinates
-  std::array<o2::math_utils::Legendre2DPolynominal, 6> mLegX;
-  std::array<o2::math_utils::Legendre2DPolynominal, 6> mLegY;
-  std::array<o2::math_utils::Legendre2DPolynominal, 6> mLegZ;
+  // Legendre polynominals to model deformations in radius; parameterized by normalized phi (u) and z (v) coordinates
+  std::array<o2::math_utils::Legendre2DPolynominal, 6> mLegendre;
 };
 
 } // namespace o2::its3::align

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

@@ -120,11 +120,6 @@ class MisAlignmentHits
 
   bool deformHit(WorkingHit::HitType t);
 
-  auto getDeformation(unsigned int id, double u, double v) const
-  {
-    return mDeformations.getDeformation(id, u, v);
-  }
-
   // Mimize function assuming a straight line
   // given in the parametric representation by y_v = t * d_x + x_s
   // assuming no offset is needed

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

@@ -16,6 +16,7 @@
 #define ITS3_MISALIGNMENTPARAMETERS_H_
 
 #include "ITS3Base/SpecsV2.h"
+#include "DetectorsCommonDataFormats/AlignParam.h"
 
 #include "TNamed.h"
 #include "TFile.h"
@@ -31,37 +32,23 @@ class MisalignmentParameters : public TNamed
 {
  public:
   MisalignmentParameters();
+  MisalignmentParameters(const MisalignmentParameters&) = default;
+  MisalignmentParameters(MisalignmentParameters&&) = delete;
+  MisalignmentParameters& operator=(MisalignmentParameters&&) = delete;
+  MisalignmentParameters& operator=(const MisalignmentParameters& o);
 
   // IO
   bool store(const std::string& file) const;
   static MisalignmentParameters* load(const std::string& file);
 
   /// Global getters
-  double getGloTransX(unsigned int detID) const { return mGloTransX[detID]; }
-  double getGloTransY(unsigned int detID) const { return mGloTransY[detID]; }
-  double getGloTransZ(unsigned int detID) const { return mGloTransZ[detID]; }
-  double getGloRotX(unsigned int detID) const { return mGloRotX[detID]; }
-  double getGloRotY(unsigned int detID) const { return mGloRotY[detID]; }
-  double getGloRotZ(unsigned int detID) const { return mGloRotZ[detID]; }
-  /// Global setters
-  void setGloTransX(unsigned int detID, double v) { mGloTransX[detID] = v; }
-  void setGloTransY(unsigned int detID, double v) { mGloTransY[detID] = v; }
-  void setGloTransZ(unsigned int detID, double v) { mGloTransZ[detID] = v; }
-  void setGloRotX(unsigned int detID, double v) { mGloRotX[detID] = v; }
-  void setGloRotY(unsigned int detID, double v) { mGloRotY[detID] = v; }
-  void setGloRotZ(unsigned int detID, double v) { mGloRotZ[detID] = v; }
-
-  /// Legendre Coeff. getters
-  const TMatrixD& getLegendreCoeffX(unsigned int sensorID) const { return mLegCoeffX[sensorID]; }
-  const TMatrixD& getLegendreCoeffY(unsigned int sensorID) const { return mLegCoeffY[sensorID]; }
-  const TMatrixD& getLegendreCoeffZ(unsigned int sensorID) const { return mLegCoeffZ[sensorID]; }
-  /// Legendre Coeff. setters
-  void setLegendreCoeffX(unsigned int sensorID, const TMatrixD& m) { setMatrix(mLegCoeffX[sensorID], m); }
-  void setLegendreCoeffY(unsigned int sensorID, const TMatrixD& m) { setMatrix(mLegCoeffY[sensorID], m); }
-  void setLegendreCoeffZ(unsigned int sensorID, const TMatrixD& m) { setMatrix(mLegCoeffZ[sensorID], m); }
+  const o2::detectors::AlignParam& getGlobal(unsigned int detID) const { return mGlobal[detID]; }
+  void setGlobal(unsigned int detID, const o2::detectors::AlignParam& param) { mGlobal[detID] = param; }
+
+  const TMatrixD& getLegendreCoeff(unsigned int sensorID) const { return mLegCoeff[sensorID]; }
+  void setLegendreCoeff(unsigned int sensorID, const TMatrixD& m) { setMatrix(mLegCoeff[sensorID], m); }
 
   void printParams(unsigned int detID) const;
-  void printLegendreParams(unsigned int sensorID) const;
 
  private:
   inline void setMatrix(TMatrixD& o, const TMatrixD& n)
@@ -71,21 +58,10 @@ class MisalignmentParameters : public TNamed
   }
 
   static constexpr unsigned int nDetectors{constants::detID::nChips}; ///! for now just the IB
+  std::array<o2::detectors::AlignParam, nDetectors> mGlobal;          ///< Array to hold the global misalignment
+  std::array<TMatrixD, constants::nSensorsIB> mLegCoeff;              // Legendre Polynominals coefficients
 
-  // Global parameters
-  std::array<double, nDetectors> mGloTransX; ///< Array to hold the global misalignment in x-direction
-  std::array<double, nDetectors> mGloTransY; ///< Array to hold the global misalignment in y-direction
-  std::array<double, nDetectors> mGloTransZ; ///< Array to hold the global misalignment in z-direction
-  std::array<double, nDetectors> mGloRotX;   ///< Array to hold the global misalignment in x-direction
-  std::array<double, nDetectors> mGloRotY;   ///< Array to hold the global misalignment in y-direction
-  std::array<double, nDetectors> mGloRotZ;   ///< Array to hold the global misalignment in z-direction
-
-  // Legendre Polynominals coefficients
-  std::array<TMatrixD, constants::nSensorsIB> mLegCoeffX;
-  std::array<TMatrixD, constants::nSensorsIB> mLegCoeffY;
-  std::array<TMatrixD, constants::nSensorsIB> mLegCoeffZ;
-
-  ClassDefOverride(MisalignmentParameters, 1);
+  ClassDefOverride(MisalignmentParameters, 2);
 };
 
 } // namespace o2::its3::align

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

@@ -27,14 +27,12 @@ void Deformations::init(const fs::path& path)
     LOGP(fatal, "File {} does not exists!", path.c_str());
   }
 
-  auto params = MisalignmentParameters::load(path.string());
+  mParams = *MisalignmentParameters::load(path.string());
   LOGP(info, "Loaded Parameters");
 
   // Set the legendre pols
   for (int iSensor{0}; iSensor < 6; ++iSensor) {
-    mLegX[iSensor] = o2::math_utils::Legendre2DPolynominal(params->getLegendreCoeffX(iSensor));
-    mLegY[iSensor] = o2::math_utils::Legendre2DPolynominal(params->getLegendreCoeffY(iSensor));
-    mLegZ[iSensor] = o2::math_utils::Legendre2DPolynominal(params->getLegendreCoeffZ(iSensor));
+    mLegendre[iSensor] = o2::math_utils::Legendre2DPolynominal(mParams.getLegendreCoeff(iSensor));
   }
 }
 

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

@@ -10,6 +10,7 @@
 // or submit itself to any jurisdiction.
 
 #include "ITS3Align/MisalignmentHits.h"
+#include "ITS3Base/SegmentationMosaix.h"
 #include "ITS3Base/ITS3Params.h"
 #include "SimConfig/DigiParams.h"
 #include "DetectorsBase/Propagator.h"
@@ -309,10 +310,8 @@ double MisAlignmentHits::StraightLine::DoEval(const double* x) const
          zline = mStart.Z() + t * mD[2];
 
   // Find the point of the deformed geometry given a certain phi' and z'
-  double xideal = mRadius * std::cos(phi), yideal = mRadius * std::sin(phi),
-         zideal = z;
-  const auto [dx, dy, dz] = mMis->getDeformation(mSensorID, nphi, nz);
-  double xdef = xideal + dx, ydef = yideal + dy, zdef = zideal + dz;
+  const double xideal = mRadius * std::cos(phi), yideal = mRadius * std::sin(phi), zideal = z;
+  const auto [xdef, ydef, zdef] = mMis->mDeformations.getDeformation(mSensorID, mRadius, phi, z, nphi, nz);
 
   // Minimize the euclidean distance of the line point and the deformed point
   return std::hypot(xline - xdef, yline - ydef, zline - zdef);
@@ -356,10 +355,7 @@ double MisAlignmentHits::Propagator::DoEval(const double* x) const
   const auto glo = trc.getXYZGlo();
 
   // Find the point of the deformed geometry given a certain phi' and z'
-  double xideal = mRadius * std::cos(phi), yideal = mRadius * std::sin(phi),
-         zideal = z;
-  const auto [dx, dy, dz] = mMis->getDeformation(mSensorID, nphi, nz);
-  double xdef = xideal + dx, ydef = yideal + dy, zdef = zideal + dz;
+  const auto [xdef, ydef, zdef] = mMis->mDeformations.getDeformation(mSensorID, mRadius, phi, z, nphi, nz);
 
   // Minimize the euclidean distance of the propagator point and the deformed point
   return std::hypot(glo.X() - xdef, glo.Y() - ydef, glo.Z() - zdef);

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

@@ -30,6 +30,20 @@ MisalignmentParameters::MisalignmentParameters()
   SetTitle("ITS3 MisalignmentParameters");
 }
 
+MisalignmentParameters& MisalignmentParameters::operator=(const MisalignmentParameters& o)
+{
+  if (this != &o) {
+    SetName(o.GetName());
+    SetTitle(o.GetTitle());
+    mGlobal = o.mGlobal;
+    for (unsigned int s{0}; s < constants::nSensorsIB; ++s) {
+      o.mLegCoeff[s].Print();
+      setMatrix(mLegCoeff[s], o.mLegCoeff[s]);
+    }
+  }
+  return *this;
+}
+
 bool MisalignmentParameters::store(const std::string& file) const
 {
   std::unique_ptr<TFile> fOut(TFile::Open(file.c_str(), "RECREATE"));
@@ -43,6 +57,7 @@ bool MisalignmentParameters::store(const std::string& file) const
 
 MisalignmentParameters* MisalignmentParameters::load(const std::string& file)
 {
+  LOGP(info, "Loading Parameters from {}", file);
   std::unique_ptr<TFile> fIn(TFile::Open(file.c_str(), "READ"));
   auto p = fIn->Get<MisalignmentParameters>("ccdb_object");
   if (p == nullptr) {
@@ -54,27 +69,13 @@ MisalignmentParameters* MisalignmentParameters::load(const std::string& file)
 void MisalignmentParameters::printParams(unsigned int detID) const
 {
   LOGP(info, "Parameters for ID={}:", detID);
-  LOGP(info, " - Global Trans: X={} Y={} Z={}", getGloTransX(detID), getGloTransY(detID), getGloTransZ(detID));
-  LOGP(info, " - Global Rots: X={} Y={} Z={}", getGloRotX(detID), getGloRotY(detID), getGloRotZ(detID));
+  const auto& glo = mGlobal[detID];
+  LOGP(info, " - Global Trans: X={} Y={} Z={}", glo.getX(), glo.getY(), glo.getZ());
+  LOGP(info, " - Global Rots: Phi={} Psi={} Theta={}", glo.getPhi(), glo.getPsi(), glo.getTheta());
   if (constants::detID::isDetITS3(detID)) {
     auto sensorID = constants::detID::getSensorID(detID);
-    LOGP(info, " - Legendre Pol X:");
-    getLegendreCoeffX(sensorID).Print();
-    LOGP(info, " - Legendre Pol Y:");
-    getLegendreCoeffY(sensorID).Print();
-    LOGP(info, " - Legendre Pol Z:");
-    getLegendreCoeffZ(sensorID).Print();
+    mLegCoeff[sensorID].Print();
   }
 }
 
-void MisalignmentParameters::printLegendreParams(unsigned int sensorID) const
-{
-  LOGP(info, " - Legendre Pol X:");
-  getLegendreCoeffX(sensorID).Print();
-  LOGP(info, " - Legendre Pol Y:");
-  getLegendreCoeffY(sensorID).Print();
-  LOGP(info, " - Legendre Pol Z:");
-  getLegendreCoeffZ(sensorID).Print();
-}
-
 } // namespace o2::its3::align

Detectors/Upgrades/ITS3/macros/align/CMakeLists.txt

@@ -15,3 +15,4 @@ its3_add_macro(TestLegendrePol.C)
 its3_add_macro(CreateMisalignmentITS3.C)
 its3_add_macro(ShowCoefficients.C)
 its3_add_macro(CheckResidualsITS3.C)
+its3_add_macro(ITS3Misaligner.C)

Detectors/Upgrades/ITS3/macros/align/CheckResidualsITS3.C

@@ -36,8 +36,12 @@
 
 #include <optional>
 #include <tuple>
+#include <filesystem>
+#include <cstdlib>
 #endif
 
+namespace fs = std::filesystem;
+
 // Refit ITS3 tracks by using the clusters from the OB only then propagate the
 // tracks to the IB clusters and calculate their residuals.
 
@@ -66,11 +70,16 @@ std::optional<Cluster> propagateTo(Track& trk, const o2::itsmft::CompClusterExt&
   ++cTotal;
   auto chipID = clus.getSensorID();
   float sigmaY2{0}, sigmaZ2{0}, sigmaYZ{0};
-  auto isITS3 = o2::its3::constants::detID::isDetITS3(chipID);
   const float alpha = o2::its::GeometryTGeo::Instance()->getSensorRefAlpha(clus.getSensorID());   // alpha for the tracking frame
   const auto locC = o2::its3::ioutils::extractClusterData(clus, pattIt, mDict, sigmaY2, sigmaZ2); // get cluster in sensor local frame with errors
-  Point3D trkC = o2::its::GeometryTGeo::Instance()->getMatrixT2L(chipID) ^ (locC);                // cluster position in the tracking frame
-  const auto gloC = o2::its::GeometryTGeo::Instance()->getMatrixL2G(chipID)(locC);                // global cluster position
+  Point3D trkC;
+  auto isITS3 = o2::its3::constants::detID::isDetITS3(chipID);
+  if (isITS3) {
+    trkC = o2::its::GeometryTGeo::Instance()->getT2LMatrixITS3(chipID, alpha) ^ (locC); // cluster position in the tracking frame
+  } else {
+    trkC = o2::its::GeometryTGeo::Instance()->getMatrixT2L(chipID) ^ (locC); // cluster position in the tracking frame
+  }
+  const auto gloC = o2::its::GeometryTGeo::Instance()->getMatrixL2G(chipID)(locC); // global cluster position
   const auto bz = o2::base::Propagator::Instance()->getNominalBz();
 
   // rotate the parameters to the tracking frame then propagate to the clusters'x
@@ -111,7 +120,7 @@ std::optional<Cluster> propagateTo(Track& trk, const o2::itsmft::CompClusterExt&
 
 void CheckResidualsITS3(bool plotOnly = false,
                         bool useITS = false,
-                        const std::string& dictFileName = "../ccdb/IT3/Calib/ClusterDictionary/snapshot.root",
+                        std::string dictFileName = "IT3dictionary.root",
                         const std::string& collisioncontextFileName = "collisioncontext.root",
                         const std::string& itsTracksFileName = "o2trac_its.root",
                         const std::string& itsClustersFileName = "o2clus_its.root",
@@ -122,6 +131,17 @@ void CheckResidualsITS3(bool plotOnly = false,
   std::array<TH3F*, 6> hIBDx, hIBDy, hIBDz, hIBDr;
 
   if (!plotOnly) {
+    if (fs::exists(dictFileName)) {
+      LOGP(info, "Using local dictionary file");
+    } else {
+      const std::string ccdb = std::getenv("ALICEO2_CCDB_LOCALCACHE");
+      dictFileName = ccdb + "/IT3/Calib/ClusterDictionary/snapshot.root";
+      if (fs::exists(dictFileName)) {
+        LOGP(info, "Loading from local ccdb cache");
+      } else {
+        LOGP(fatal, "Provide dictionary!");
+      }
+    }
     mDict = new o2::its3::TopologyDictionary(dictFileName);
 
     o2::base::GeometryManager::loadGeometry(geomFileName);