Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:33 

0001 /****************************************************************************
0002 * Authors: 
0003 *  Jan Kašpar (jan.kaspar@gmail.com) 
0004 ****************************************************************************/
0005 
0006 #include "CalibPPS/AlignmentRelative/interface/IdealResult.h"
0007 
0008 #include "CalibPPS/AlignmentRelative/interface/AlignmentTask.h"
0009 
0010 #include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
0011 #include "Geometry/Records/interface/VeryForwardMisalignedGeometryRecord.h"
0012 
0013 #include "TMatrixD.h"
0014 #include "TVectorD.h"
0015 
0016 using namespace std;
0017 
0018 //----------------------------------------------------------------------------------------------------
0019 
0020 IdealResult::IdealResult(const edm::ParameterSet &gps, AlignmentTask *_t) : AlignmentAlgorithm(gps, _t) {}
0021 
0022 //----------------------------------------------------------------------------------------------------
0023 
0024 void IdealResult::begin(const CTPPSGeometry *geometryReal, const CTPPSGeometry *geometryMisaligned) {
0025   gReal = geometryReal;
0026   gMisaligned = geometryMisaligned;
0027 }
0028 //----------------------------------------------------------------------------------------------------
0029 
0030 unsigned int IdealResult::solve(const std::vector<AlignmentConstraint> &constraints,
0031                                 std::map<unsigned int, AlignmentResult> &results,
0032                                 TDirectory *dir) {
0033   /* 
0034   STRATEGY:
0035 
0036   1) Determine true misalignment as misalign_geometry - real_geometry.
0037 
0038   2) Represent the true misalignment as the "chi" vector (cf. Jan algorithm), denoted chi_tr.
0039 
0040   3) Find the expected result of track-based alignment, subject to the given constraints. Formally this corresponds to finding
0041        chi_exp = chi_tr + I * Lambda
0042      where the I matrix contains the "inaccessible alignment modes" as columns and Lambda is a vector of parameters. The constraints are given by
0043        C^T * chi_exp = V
0044      Since the problem may be generally overconstrained, it is better to formulate it as search for Lambda which minimises
0045        |C^ * chi_exp - V|^2
0046      This gives
0047        Lambda = (A^T * A)^-1 * A^T * b,  A = C^T * I,  b = V - C^T * chi_tr
0048   */
0049 
0050   results.clear();
0051 
0052   // determine dimension and offsets
0053   unsigned int dim = 0;
0054   vector<unsigned int> offsets;
0055   map<unsigned int, unsigned int> offset_map;
0056   for (unsigned int qci = 0; qci < task->quantityClasses.size(); qci++) {
0057     offsets.push_back(dim);
0058     offset_map[task->quantityClasses[qci]] = dim;
0059     dim += task->quantitiesOfClass(task->quantityClasses[qci]);
0060   }
0061 
0062   // collect true misalignments
0063   TVectorD chi_tr(dim);
0064 
0065   for (const auto &dit : task->geometry.getSensorMap()) {
0066     unsigned int detId = dit.first;
0067 
0068     const DetGeomDesc *real = gReal->sensor(detId);
0069     const DetGeomDesc *misal = gMisaligned->sensor(detId);
0070 
0071     // extract shift
0072     const auto shift = misal->translation() - real->translation();
0073 
0074     // extract rotation around z
0075     const auto rotation = misal->rotation() * real->rotation().Inverse();
0076 
0077     double r_xx, r_xy, r_xz;
0078     double r_yx, r_yy, r_yz;
0079     double r_zx, r_zy, r_zz;
0080     rotation.GetComponents(r_xx, r_xy, r_xz, r_yx, r_yy, r_yz, r_zx, r_zy, r_zz);
0081 
0082     if (std::abs(r_zz - 1.) > 1E-5)
0083       throw cms::Exception("PPS") << "IdealResult::Solve: only rotations about z are supported.";
0084 
0085     double rot_z = atan2(r_yx, r_xx);
0086 
0087     const auto &geom = task->geometry.get(detId);
0088 
0089     for (unsigned int qci = 0; qci < task->quantityClasses.size(); ++qci) {
0090       const auto &qc = task->quantityClasses[qci];
0091       signed int idx = task->getQuantityIndex(qc, detId);
0092       if (idx < 0)
0093         continue;
0094 
0095       double v = 0.;
0096 
0097       if (qc == AlignmentTask::qcShR1) {
0098         const auto &d = geom.getDirectionData(1);
0099         v = shift.x() * d.dx + shift.y() * d.dy + shift.z() * d.dz;
0100       }
0101 
0102       if (qc == AlignmentTask::qcShR2) {
0103         const auto &d = geom.getDirectionData(2);
0104         v = shift.x() * d.dx + shift.y() * d.dy + shift.z() * d.dz;
0105       }
0106 
0107       if (qc == AlignmentTask::qcRotZ)
0108         v = rot_z;
0109 
0110       chi_tr(offsets[qci] + idx) = v;
0111     }
0112   }
0113 
0114   // build list of "inaccessible" modes
0115   vector<TVectorD> inaccessibleModes;
0116 
0117   if (task->resolveShR) {
0118     TVectorD fm_ShX_gl(dim);
0119     fm_ShX_gl.Zero();
0120     TVectorD fm_ShX_lp(dim);
0121     fm_ShX_lp.Zero();
0122     TVectorD fm_ShY_gl(dim);
0123     fm_ShY_gl.Zero();
0124     TVectorD fm_ShY_lp(dim);
0125     fm_ShY_lp.Zero();
0126 
0127     for (const auto &sp : task->geometry.getSensorMap()) {
0128       CTPPSDetId senId(sp.first);
0129       const auto &geom = sp.second;
0130 
0131       if (senId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0132         signed int qIndex = task->getQuantityIndex(AlignmentTask::qcShR2, senId);
0133 
0134         const double d2x = geom.getDirectionData(2).dx;
0135         const double d2y = geom.getDirectionData(2).dy;
0136 
0137         const auto &offset2 = offset_map[AlignmentTask::qcShR2];
0138         fm_ShX_gl(offset2 + qIndex) = d2x;
0139         fm_ShX_lp(offset2 + qIndex) = d2x * geom.z;
0140         fm_ShY_gl(offset2 + qIndex) = d2y;
0141         fm_ShY_lp(offset2 + qIndex) = d2y * geom.z;
0142       }
0143 
0144       if (senId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0145         const signed int qIndex1 = task->getQuantityIndex(AlignmentTask::qcShR1, senId);
0146         const signed int qIndex2 = task->getQuantityIndex(AlignmentTask::qcShR2, senId);
0147 
0148         const double d1x = geom.getDirectionData(1).dx;
0149         const double d1y = geom.getDirectionData(1).dy;
0150         const double d2x = geom.getDirectionData(2).dx;
0151         const double d2y = geom.getDirectionData(2).dy;
0152 
0153         const auto &offset1 = offset_map[AlignmentTask::qcShR1];
0154         fm_ShX_gl(offset1 + qIndex1) = d1x;
0155         fm_ShX_lp(offset1 + qIndex1) = d1x * geom.z;
0156         fm_ShY_gl(offset1 + qIndex1) = d1y;
0157         fm_ShY_lp(offset1 + qIndex1) = d1y * geom.z;
0158 
0159         const auto &offset2 = offset_map[AlignmentTask::qcShR2];
0160         fm_ShX_gl(offset2 + qIndex2) = d2x;
0161         fm_ShX_lp(offset2 + qIndex2) = d2x * geom.z;
0162         fm_ShY_gl(offset2 + qIndex2) = d2y;
0163         fm_ShY_lp(offset2 + qIndex2) = d2y * geom.z;
0164       }
0165     }
0166 
0167     inaccessibleModes.push_back(fm_ShX_gl);
0168     inaccessibleModes.push_back(fm_ShX_lp);
0169     inaccessibleModes.push_back(fm_ShY_gl);
0170     inaccessibleModes.push_back(fm_ShY_lp);
0171   }
0172 
0173   if (task->resolveRotZ) {
0174     TVectorD fm_RotZ_gl(dim);
0175     fm_RotZ_gl.Zero();
0176     TVectorD fm_RotZ_lp(dim);
0177     fm_RotZ_lp.Zero();
0178 
0179     for (const auto &sp : task->geometry.getSensorMap()) {
0180       CTPPSDetId senId(sp.first);
0181       const auto &geom = sp.second;
0182 
0183       for (int m = 0; m < 2; ++m) {
0184         double rho = 0.;
0185         TVectorD *fmp = nullptr;
0186         if (m == 0) {
0187           rho = 1.;
0188           fmp = &fm_RotZ_gl;
0189         }
0190         if (m == 1) {
0191           rho = geom.z;
0192           fmp = &fm_RotZ_lp;
0193         }
0194         TVectorD &fm = *fmp;
0195 
0196         const signed int qIndex = task->getQuantityIndex(AlignmentTask::qcRotZ, senId);
0197         const auto &offset = offset_map[AlignmentTask::qcRotZ];
0198         fm(offset + qIndex) = rho;
0199 
0200         const double as_x = -rho * geom.sy;
0201         const double as_y = +rho * geom.sx;
0202 
0203         if (senId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0204           const double d2x = geom.getDirectionData(2).dx;
0205           const double d2y = geom.getDirectionData(2).dy;
0206 
0207           const signed int qIndex2 = task->getQuantityIndex(AlignmentTask::qcShR2, senId);
0208           const auto &offset2 = offset_map[AlignmentTask::qcShR2];
0209           fm(offset2 + qIndex2) = d2x * as_x + d2y * as_y;
0210         }
0211 
0212         if (senId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0213           const double d1x = geom.getDirectionData(1).dx;
0214           const double d1y = geom.getDirectionData(1).dy;
0215           const double d2x = geom.getDirectionData(2).dx;
0216           const double d2y = geom.getDirectionData(2).dy;
0217 
0218           const signed int qIndex1 = task->getQuantityIndex(AlignmentTask::qcShR1, senId);
0219           const auto &offset1 = offset_map[AlignmentTask::qcShR1];
0220           fm(offset1 + qIndex1) = d1x * as_x + d1y * as_y;
0221 
0222           const signed int qIndex2 = task->getQuantityIndex(AlignmentTask::qcShR2, senId);
0223           const auto &offset2 = offset_map[AlignmentTask::qcShR2];
0224           fm(offset2 + qIndex2) = d2x * as_x + d2y * as_y;
0225         }
0226       }
0227     }
0228 
0229     inaccessibleModes.push_back(fm_RotZ_gl);
0230     inaccessibleModes.push_back(fm_RotZ_lp);
0231   }
0232 
0233   // build matrices and vectors
0234   TMatrixD C(dim, constraints.size());
0235   TVectorD V(constraints.size());
0236   for (unsigned int i = 0; i < constraints.size(); i++) {
0237     V(i) = constraints[i].val;
0238 
0239     unsigned int offset = 0;
0240     for (auto &quantityClass : task->quantityClasses) {
0241       const TVectorD &cv = constraints[i].coef.find(quantityClass)->second;
0242       for (int k = 0; k < cv.GetNrows(); k++) {
0243         C[offset][i] = cv[k];
0244         offset++;
0245       }
0246     }
0247   }
0248 
0249   TMatrixD I(dim, inaccessibleModes.size());
0250   for (unsigned int i = 0; i < inaccessibleModes.size(); ++i) {
0251     for (int j = 0; j < inaccessibleModes[i].GetNrows(); ++j)
0252       I(j, i) = inaccessibleModes[i](j);
0253   }
0254 
0255   // determine expected track-based alignment result
0256   TMatrixD CT(TMatrixD::kTransposed, C);
0257   TMatrixD CTI(CT * I);
0258 
0259   const TMatrixD &A = CTI;
0260   TMatrixD AT(TMatrixD::kTransposed, A);
0261   TMatrixD ATA(AT * A);
0262   TMatrixD ATA_inv(TMatrixD::kInverted, ATA);
0263 
0264   TVectorD b = V - CT * chi_tr;
0265 
0266   TVectorD La(ATA_inv * AT * b);
0267 
0268   TVectorD chi_exp(chi_tr + I * La);
0269 
0270   // save result
0271   for (const auto &dit : task->geometry.getSensorMap()) {
0272     AlignmentResult r;
0273 
0274     for (unsigned int qci = 0; qci < task->quantityClasses.size(); ++qci) {
0275       const auto &qc = task->quantityClasses[qci];
0276       const auto idx = task->getQuantityIndex(qc, dit.first);
0277       if (idx < 0)
0278         continue;
0279 
0280       const auto &v = chi_exp(offsets[qci] + idx);
0281 
0282       switch (qc) {
0283         case AlignmentTask::qcShR1:
0284           r.setShR1(v);
0285           break;
0286         case AlignmentTask::qcShR2:
0287           r.setShR2(v);
0288           break;
0289         case AlignmentTask::qcShZ:
0290           r.setShZ(v);
0291           break;
0292         case AlignmentTask::qcRotZ:
0293           r.setRotZ(v);
0294           break;
0295       }
0296     }
0297 
0298     results[dit.first] = r;
0299   }
0300 
0301   return 0;
0302 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!