Finally reasonable results.

cms-sw · Sep 4, 2020 · fd58ad1 · fd58ad1
1 parent 1266951
commit fd58ad1
Showing 1 changed file with 39 additions and 59 deletions.
diff --git a/Alignment/PPSTrackBased/src/IdealResult.cc b/Alignment/PPSTrackBased/src/IdealResult.cc
@@ -42,6 +42,23 @@ vector<SingularMode> IdealResult::analyze()
 unsigned int IdealResult::solve(const std::vector<AlignmentConstraint> &constraints,
   std::map<unsigned int, AlignmentResult> &results, TDirectory *dir)
 {
+  /* 
+  STRATEGY:
+
+  1) Determine true misalignment as misalign_geometry - real_geometry.
+
+  2) Represent the true misalignment as the "chi" vector (cf. Jan algorithm), denoted chi_tr.
+
+  3) Find the expected result of track-based alignment, subject to the given constraints. Formally this corresponds to finding
+       chi_exp = chi_tr + I * Lambda
+     where the I matrix contains the "inaccessible alignment modes" as columns and Lambda is a vector of parameters. The constraints are given by
+       C^T * chi_exp = V
+     Since the problem may be generally overconstrained, it is better to formulate it as search for Lambda which minimises
+       |C^ * chi_exp - V|^2
+     This gives
+       Lambda = (A^T * A)^-1 * A^T * b,  A = C^T * I,  b = V - C^T * chi_tr
+  */
+
   results.clear();
 
   // determine dimension and offsets
@@ -227,49 +244,13 @@ unsigned int IdealResult::solve(const std::vector<AlignmentConstraint> &constrai
     inaccessibleModes.push_back(fm_RotZ_lp);
   }
 
-  // ortonormalise inaccessible modes (Gramm-Schmidt)
-  vector<TVectorD> inaccessibleModesON;
-  for (const auto &fm : inaccessibleModes)
-  {
-    unsigned int i = inaccessibleModesON.size();
-
-    auto q = fm;
-
-    for (unsigned int k = 0; k < i; ++k)
-    {
-      double d = 0;
-      for (unsigned int idx = 0; idx < dim; ++idx)
-        d += inaccessibleModes[i][idx] * inaccessibleModesON[k][idx];
-
-      for (unsigned int idx = 0; idx < dim; ++idx)
-        q(idx) -= d * inaccessibleModesON[k][idx];
-    }
-
-    q *= 1./sqrt(q.Norm2Sqr());
-    inaccessibleModesON.push_back(q);
-  }
-
-  // simulate S matrix as projector to accessible modes
-  TMatrixDSym S(dim);
-  S.Zero();
-  for (unsigned int i = 0; i < dim; i++)
-      S(i, i) = 1.;
-  for (const auto &fm : inaccessibleModesON)
-  {
-    for (unsigned int i = 0; i < dim; i++)
-    {
-      for (unsigned int j = 0; j < dim; j++)
-        S(i, j) -= fm(i) * fm(j);
-    }
-  }
-
-  // simulate the M vector
-  TVectorD M(S * chi_tr);
-
-  // build constraint matrix
+  // build matrices and vectors
   TMatrixD C(dim, constraints.size());
+  TVectorD V(constraints.size());
   for (unsigned int i = 0; i < constraints.size(); i++)
   {
+    V(i) = constraints[i].val;
+
     unsigned int offset = 0;
     for (auto &quantityClass : task->quantityClasses)
     {
@@ -282,28 +263,27 @@ unsigned int IdealResult::solve(const std::vector<AlignmentConstraint> &constrai
     }
   }
 
-  // simulate CS matrix
-  TMatrixDSym CS(dim + constraints.size());
-  CS.Zero();
-  for (unsigned int i = 0; i < dim; i++)
-      CS(i, i) = S(i, i);
-
-  for (unsigned int i = 0; i < constraints.size(); i++)
+  TMatrixD I(dim, inaccessibleModes.size());
+  for (unsigned int i = 0; i < inaccessibleModes.size(); ++i)
   {
-    for (unsigned int j = 0; j < dim; j++)
-      CS[j][dim + i] = CS[dim + i][j] =  C(j, i);
+    for (int j = 0; j < inaccessibleModes[i].GetNrows(); ++j)
+      I(j, i) = inaccessibleModes[i](j);
   }
 
-  // build MV vector
-  TVectorD MV(dim + constraints.size());
-  for (unsigned int i = 0; i < dim; i++)
-    MV(i) = M(i);
-  for (unsigned int i = 0; i < constraints.size(); i++)
-    MV[dim + i] = constraints[i].val;
+  // determine expected track-based alignment result
+  TMatrixD CT(TMatrixD::kTransposed, C);
+  TMatrixD CTI(CT * I);
+
+  TMatrixD A = CTI;
+  TMatrixD AT(TMatrixD::kTransposed, A);
+  TMatrixD ATA(AT*A);
+  TMatrixD ATA_inv(TMatrixD::kInverted, ATA);
+
+  TVectorD b = V - CT * chi_tr;
+
+  TVectorD La(ATA_inv * AT * b);
 
-  // solve: CS * R = MV
-  TMatrixD CS_inv(TMatrixD::kInverted, CS);
-  TVectorD R(CS_inv * MV);
+  TVectorD chi_exp(chi_tr + I * La);
 
   // save result
   for (const auto &dit : task->geometry.getSensorMap())
@@ -317,7 +297,7 @@ unsigned int IdealResult::solve(const std::vector<AlignmentConstraint> &constrai
       if (idx < 0)
         continue;
 
-      const auto &v = R(offsets[qci] + idx);
+      const auto &v = chi_exp(offsets[qci] + idx);
 
       switch (qc)
       {