align

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

Author: f3sch

Detectors/Upgrades/ITS3/alignment/CMakeLists.txt

@@ -27,6 +27,7 @@ o2_add_library(ITS3Align
                                      O2::GlobalTrackingWorkflowHelpers
                                      O2::DataFormatsGlobalTracking
                                      O2::DetectorsVertexing
+                                     nlohmann_json::nlohmann_json
                                      GBL::GBL)
 if (OpenMP_CXX_FOUND)
     target_compile_definitions(${targetName} PRIVATE WITH_OPENMP)

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

@@ -36,13 +36,17 @@ struct AlignmentParams : public o2::conf::ConfigurableParamHelper<AlignmentParam
   float extraClsErrY[7] = {0};
   float extraClsErrZ[7] = {0};
 
+  // misalignment
+  bool doMisalignment = false;   // simulate Legendre deformation on ITS3 layers
+  bool doMisalignmentRB = false; // simulate rigid body misalignment on ITS3 layers
+  std::string misAlgJson;        // JSON file with deformation and/or rigid body params
+
   // 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.2; // max tolerance of diagonal elements away from 1
 
   // MillePede output
   std::string milleBinFile = "mp2data.bin";

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

@@ -23,6 +23,7 @@
 #include <GblPoint.h>
 #include <GblMeasurement.h>
 #include <MilleBinary.h>
+#include <nlohmann/json.hpp>
 
 #include "Framework/ConfigParamRegistry.h"
 #include "Framework/DataProcessorSpec.h"
@@ -47,6 +48,7 @@
 #include "ITS3Align/AlignmentTypes.h"
 #include "ITS3Align/AlignmentHierarchy.h"
 #include "ITS3Align/AlignmentSensors.h"
+#include "MathUtils/LegendrePols.h"
 
 namespace o2::its3::align
 {
@@ -97,6 +99,11 @@ class AlignmentSpec final : public Task
   // build track to vertex association
   void buildT2V();
 
+  // apply some misalignment on inner ITS3 layers
+  // it can happen that a measurement is pushed outside of
+  // ITS3 acceptance so false is to discard track
+  bool applyMisalignment(Eigen::Vector2d& res, const FrameInfoExt& frame, const TrackD& wTrk, size_t iTrk);
+
   // Steering debug
   enum class OutputOpt : uint8_t {
     VerboseGBL = 0,
@@ -119,10 +126,11 @@ class AlignmentSpec final : public Task
   AlignableVolume::SensorMapping mChip2Hiearchy; // global label mapping to leaves in the tree
   bool mUseMC{false};
   bool mWithPV{false};
-  FrameInfoExt mPVInfo;
   GTrackID::mask_t mTracksSrc;
   int mNThreads{1};
   const AlignmentParams* mParams{nullptr};
+  std::array<o2::math_utils::Legendre2DPolynominal, 6> mDeformations;       // one per sensorID (0-5)
+  std::array<Eigen::Matrix<double, 6, 1>, 6> mRigidBodyParams;             // (dx,dy,dz,rx,ry,rz) in LOC per sensorID
 };
 
 void AlignmentSpec::init(InitContext& ic)
@@ -264,8 +272,9 @@ void AlignmentSpec::process()
         res << frame.positionTrackingFrame[0] - wTrk.getY(), frame.positionTrackingFrame[1] - wTrk.getZ();
 
         // here we can apply some misalignment on the measurment
-        if (constants::detID::getSensorID(frame.sens) == 2) {
-          res[0] += 0.01;
+        if (!applyMisalignment(res, frame, wTrk, iTrk)) {
+          failed = true;
+          break;
         }
 
         prec << 1. / resTrack.points[ip].sig2y, 1. / resTrack.points[ip].sig2z;
@@ -274,7 +283,7 @@ void AlignmentSpec::process()
         if (msErr > mParams->minMS && ip < np - 1) {
           Eigen::Vector2d scat(0., 0.), scatPrec = Eigen::Vector2d::Constant(1. / (msErr * msErr));
           point.addScatterer(scat, scatPrec);
-          lt.clearFast(); // only clear if accounted
+          lt.clearFast(); // clear if accounted
         }
 
         if (frame.lr >= 0) {
@@ -356,9 +365,8 @@ void AlignmentSpec::process()
               LOGP(info, "GBL FIT chi2 {} ndf {}", chi2, ndf);
               traj.printTrajectory(5);
             }
-            if (chi2 / ndf > mParams->maxChi2Ndf) {
+            if (chi2 / ndf > mParams->maxChi2Ndf && cGBLChi2Rej < 10) {
               LOGP(error, "GBL fit exceeded red chi2 {}", chi2 / ndf);
-              ++cGBLChi2Rej;
               if (std::abs(resTrack.kfFit.chi2Ndf - 1) < 0.02) {
                 LOGP(error, "\tGBL is far away from good KF fit!!!!");
                 continue;
@@ -431,6 +439,39 @@ void AlignmentSpec::updateTimeDependentParams(ProcessingContext& pc)
     mParams->printKeyValues(true, true);
 
     buildHierarchy();
+
+    if (mParams->doMisalignment || 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->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<int>();
+          if (mParams->doMisalignment && item.contains("matrix")) {
+            auto v = item["matrix"].get<std::vector<std::vector<double>>>();
+            TMatrixD m(v.size(), v[0].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);
+          }
+          if (mParams->doMisalignmentRB && item.contains("rigidBody")) {
+            auto rb = item["rigidBody"].get<std::vector<double>>();
+            for (int k = 0; k < 6 && k < (int)rb.size(); ++k) {
+              mRigidBodyParams[id](k) = rb[k];
+            }
+          }
+        }
+      }
+    }
   }
 }
 
@@ -536,12 +577,6 @@ bool AlignmentSpec::getTransportJacobian(const TrackD& track, double xTo, double
     LOGP(error, "Near jacobian idendity for taking track from {} to {}", track.getX(), xTo);
     return false;
   }
-  // only check for elements where there is no change expected q/pt, tgl
-  if (std::abs(1. - jac(4, 4)) > mParams->ridderMaxJacDiagTol || std::abs(1. - jac(3, 3)) > mParams->ridderMaxJacDiagTol) {
-    LOGP(error, "Diagonal element not expected to change in propagation jacobian is to far away from 1");
-    LOG(debug) << jac;
-    return false;
-  }
 
   return true;
 }
@@ -584,6 +619,7 @@ bool AlignmentSpec::prepareITSTrack(int iTrk, const o2::its::TrackITS& itsTrack,
     return 1;
   };
 
+  FrameInfoExt pvInfo;
   if (mWithPV) { // add PV as constraint
     const int iPV = mT2PV[iTrk];
     if (iPV < 0) {
@@ -594,18 +630,18 @@ bool AlignmentSpec::prepareITSTrack(int iTrk, const o2::its::TrackITS& itsTrack,
     if (!prop->propagateToDCA(pv, tmp, bz)) {
       return false;
     }
-    mPVInfo.alpha = (float)tmp.getAlpha();
+    pvInfo.alpha = (float)tmp.getAlpha();
     double ca{0}, sa{0};
-    o2::math_utils::bringTo02Pid(mPVInfo.alpha);
-    o2::math_utils::sincosd(mPVInfo.alpha, sa, ca);
-    mPVInfo.x = tmp.getX();
-    mPVInfo.positionTrackingFrame[0] = -pv.getX() * sa + pv.getY() * ca;
-    mPVInfo.positionTrackingFrame[1] = pv.getZ();
-    mPVInfo.covarianceTrackingFrame[0] = 0.5 * (pv.getSigmaX2() + pv.getSigmaY2());
-    mPVInfo.covarianceTrackingFrame[2] = pv.getSigmaY2();
-    mPVInfo.sens = -1;
-    mPVInfo.lr = -1;
-    frameArr[0] = &mPVInfo;
+    o2::math_utils::bringToPMPid(pvInfo.alpha);
+    o2::math_utils::sincosd(pvInfo.alpha, sa, ca);
+    pvInfo.x = tmp.getX();
+    pvInfo.positionTrackingFrame[0] = -pv.getX() * sa + pv.getY() * ca;
+    pvInfo.positionTrackingFrame[1] = pv.getZ();
+    pvInfo.covarianceTrackingFrame[0] = 0.5 * (pv.getSigmaX2() + pv.getSigmaY2());
+    pvInfo.covarianceTrackingFrame[2] = pv.getSigmaY2();
+    pvInfo.sens = -1;
+    pvInfo.lr = -1;
+    frameArr[0] = &pvInfo;
   }
 
   // collect all track clusters to array, placing them to layer+1 slot
@@ -723,6 +759,105 @@ void AlignmentSpec::buildT2V()
   }
 }
 
+bool AlignmentSpec::applyMisalignment(Eigen::Vector2d& res, const FrameInfoExt& frame, const TrackD& wTrk, size_t iTrk)
+{
+  if (!constants::detID::isDetITS3(frame.sens)) {
+    return true;
+  }
+
+  const int sensorID = constants::detID::getSensorID(frame.sens);
+  const int layerID = constants::detID::getDetID2Layer(frame.sens);
+
+  // --- Legendre deformation (non-rigid-body) ---
+  if (mParams->doMisalignment && mIT3Dict && mUseMC) {
+    const auto prop = o2::base::PropagatorD::Instance();
+    // compute normalized (u,v) in [-1,1] from global position
+    auto computeUV = [](double gloX, double gloY, double gloZ, int sensorID, double radius) -> std::pair<double, double> {
+      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 lbl = mRecoData->getITSTracksMCLabels()[iTrk];
+    const auto mcTrk = mcReader->getTrack(lbl);
+    if (!mcTrk) {
+      return false;
+    }
+    std::array<double, 3> xyz{mcTrk->GetStartVertexCoordinatesX(), mcTrk->GetStartVertexCoordinatesY(), mcTrk->GetStartVertexCoordinatesZ()};
+    std::array<double, 3> pxyz{mcTrk->GetStartVertexMomentumX(), mcTrk->GetStartVertexMomentumY(), mcTrk->GetStartVertexMomentumZ()};
+    TParticlePDG* pPDG = TDatabasePDG::Instance()->GetParticle(mcTrk->GetPdgCode());
+    if (!pPDG) {
+      return false;
+    }
+    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.) {
+      return false;
+    }
+
+    res[0] += dy;
+    res[1] += dz;
+  }
+
+  // --- Rigid body misalignment ---
+  // Uses the same derivative chain as GBL global derivatives:
+  //   dres/da_LOC = dres/da_TRK * J_L2T
+  // so the residual shift is der_LOC * delta_a
+  if (mParams->doMisalignmentRB) {
+    GlobalLabel lbl(0, frame.sens, true);
+    if (mChip2Hiearchy.find(lbl) == mChip2Hiearchy.end()) {
+      return true; // sensor not in hierarchy, skip
+    }
+    const auto* sensorVol = mChip2Hiearchy.at(lbl);
+
+    // derivative in TRK frame (3x6: rows = dy, dz, dsnp)
+    Matrix36 der = getRigidBodyDerivatives(wTrk);
+
+    // transform TRK -> LOC
+    const double posTrk[3] = {frame.x, 0., 0.};
+    double posLoc[3];
+    sensorVol->getT2L().LocalToMaster(posTrk, posLoc);
+    Matrix66 jacL2T;
+    sensorVol->computeJacobianL2T(posLoc, jacL2T);
+    der *= jacL2T;
+
+    // apply: delta_res = der(row 0,1) * delta_a_LOC
+    Eigen::Vector3d shift = der * mRigidBodyParams[sensorID];
+    res[0] += shift[0]; // dy
+    res[1] += shift[1]; // dz
+  }
+
+  return true;
+}
+
 void AlignmentSpec::endOfStream(EndOfStreamContext& /*ec*/)
 {
   mDBGOut->Close();

Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx

@@ -1103,7 +1103,7 @@ 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 δy = dydx·h(u,v), δz = dzdx·h(u,v) — first Newton step.
+    // then compute dy = dydx*h(u,v), dz = dzdx*h(u,v) - first Newton step.
     const auto mcTrk = mMCReader.getTrack(lbl);
     if (!mcTrk) {
       continue;