AlignmentRelative/src/Utilities.cc

0001 /****************************************************************************
0002 * Authors:
0003 *  Jan Kašpar (jan.kaspar@gmail.com)
0004 ****************************************************************************/
0005
0006 #include "CalibPPS/AlignmentRelative/interface/Utilities.h"
0007
0008 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0009 #include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h"
0010 #include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h"
0011 #include "DataFormats/CTPPSDetId/interface/CTPPSPixelDetId.h"
0012
0013 #include "CalibPPS/AlignmentRelative/interface/AlignmentGeometry.h"
0014 #include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsData.h"
0015
0016 #include "TVectorD.h"
0017 #include "TMatrixD.h"
0018
0019 #include <map>
0020 #include <set>
0021
0022 using namespace std;
0023
0024 //----------------------------------------------------------------------------------------------------
0025
0026 void printId(unsigned int id) {
0027   CTPPSDetId detId(id);
0028
0029   if (detId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0030     TotemRPDetId stDetId(id);
0031     printf("strip %u (%3u.%u)", id, 100 * stDetId.arm() + 10 * stDetId.station() + stDetId.rp(), stDetId.plane());
0032   }
0033
0034   if (detId.subdetId() == CTPPSDetId::sdTimingDiamond) {
0035     CTPPSDiamondDetId diDetId(id);
0036     printf("dimnd %u (%3u.%u)", id, 100 * diDetId.arm() + 10 * diDetId.station() + diDetId.rp(), diDetId.plane());
0037   }
0038
0039   if (detId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0040     CTPPSPixelDetId piDetId(id);
0041     printf("pixel %u (%3u.%u)", id, 100 * piDetId.arm() + 10 * piDetId.station() + piDetId.rp(), piDetId.plane());
0042   }
0043 }
0044
0045 //----------------------------------------------------------------------------------------------------
0046
0047 void print(TMatrixD &m, const char *label, bool mathematicaFormat) {
0048   if (mathematicaFormat) {
0049     printf("{");
0050     for (int i = 0; i < m.GetNrows(); i++) {
0051       if (i > 0)
0052         printf(", ");
0053       printf("{");
0054       for (int j = 0; j < m.GetNcols(); j++) {
0055         if (j > 0)
0056           printf(", ");
0057         printf("%.3f", m[i][j]);
0058       }
0059       printf("}");
0060     }
0061     printf("}\n");
0062     return;
0063   }
0064
0065   if (label)
0066     printf("\n%s\n", label);
0067
0068   printf("    | ");
0069   for (int j = 0; j < m.GetNcols(); j++)
0070     printf(" %9i", j);
0071   printf("\n------");
0072   for (int j = 0; j < m.GetNcols(); j++)
0073     printf("----------");
0074   printf("\n");
0075
0076   for (int i = 0; i < m.GetNrows(); i++) {
0077     printf("%3i | ", i);
0078     for (int j = 0; j < m.GetNcols(); j++) {
0079       double v = m[i][j];
0080       if (fabs(v) >= 1E4)
0081         printf(" %+9.2E", v);
0082       else if (fabs(v) > 1E-6)
0083         printf(" %+9.2E", v);
0084       else
0085         printf("         0");
0086     }
0087     printf("\n");
0088   }
0089 }
0090
0091 //----------------------------------------------------------------------------------------------------
0092
0093 void factorRPFromSensorCorrections(const CTPPSRPAlignmentCorrectionsData &inputAlignments,
0094                                    CTPPSRPAlignmentCorrectionsData &expandedAlignments,
0095                                    CTPPSRPAlignmentCorrectionsData &factoredAlignments,
0096                                    const AlignmentGeometry &geometry,
0097                                    bool equalWeights,
0098                                    unsigned int verbosity) {
0099   // clean first
0100   expandedAlignments.clear();
0101   factoredAlignments.clear();
0102
0103   // save full sensor alignments
0104   map<unsigned int, CTPPSRPAlignmentCorrectionData> fullAlignments;
0105   map<unsigned int, set<unsigned int> > sensorsPerRP;
0106   for (auto it : inputAlignments.getSensorMap()) {
0107     const auto &sensorId = it.first;
0108
0109     // with useRPErrors=false the only the sensor uncertainties (coming from the last analysis step) will be used
0110     fullAlignments[sensorId] = inputAlignments.getFullSensorCorrection(sensorId, false);
0111
0112     sensorsPerRP[CTPPSDetId(sensorId).rpId()].insert(sensorId);
0113   }
0114
0115   // convert full alignments to expandedAlignments
0116   for (const auto &it : fullAlignments) {
0117     expandedAlignments.setSensorCorrection(it.first, it.second);
0118   }
0119
0120   // do the factorization RP per RP
0121   for (const auto &rpit : sensorsPerRP) {
0122     CTPPSDetId rpId(rpit.first);
0123     const set<unsigned int> &sensors = rpit.second;
0124
0125     if (verbosity)
0126       printf("* processing RP %u (%u)\n", rpit.first, 100 * rpId.arm() + 10 * rpId.station() + rpId.rp());
0127
0128     // determine number of constraints
0129     unsigned int n_constraints = 0;
0130     for (const auto &senId : sensors) {
0131       CTPPSDetId detId(senId);
0132
0133       if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip)
0134         n_constraints += 1;
0135
0136       if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel)
0137         n_constraints += 2;
0138     }
0139
0140     // build matrices
0141     TMatrixD B(n_constraints, 2), Vi(n_constraints, n_constraints), VarM(n_constraints, n_constraints);
0142     TVectorD M(n_constraints);
0143
0144     double sw2_sh_z = 0., svw2_sh_z = 0., su2w4_sh_z = 0.;
0145     double sw2_rot_x = 0., svw2_rot_x = 0., su2w4_rot_x = 0.;
0146     double sw2_rot_y = 0., svw2_rot_y = 0., su2w4_rot_y = 0.;
0147     double sw2_rot_z = 0., svw2_rot_z = 0., su2w4_rot_z = 0.;
0148
0149     unsigned int idx = 0;
0150     for (const auto &senId : sensors) {
0151       CTPPSDetId detId(senId);
0152
0153       const double v_sh_z = fullAlignments[senId].getShZ();
0154       const double u_sh_z = (fullAlignments[senId].getShZUnc() > 0.) ? fullAlignments[senId].getShZUnc() : 1.;
0155       const double w_sh_z = (equalWeights) ? 1. : 1. / u_sh_z;
0156       sw2_sh_z += w_sh_z * w_sh_z;
0157       svw2_sh_z += v_sh_z * w_sh_z * w_sh_z;
0158       su2w4_sh_z += u_sh_z * u_sh_z * w_sh_z * w_sh_z * w_sh_z * w_sh_z;
0159
0160       const double v_rot_x = fullAlignments[senId].getRotX();
0161       const double u_rot_x = (fullAlignments[senId].getRotXUnc() > 0.) ? fullAlignments[senId].getRotXUnc() : 1.;
0162       const double w_rot_x = (equalWeights) ? 1. : 1. / u_rot_x;
0163       sw2_rot_x += w_rot_x * w_rot_x;
0164       svw2_rot_x += v_rot_x * w_rot_x * w_rot_x;
0165       su2w4_rot_x += u_rot_x * u_rot_x * w_rot_x * w_rot_x * w_rot_x * w_rot_x;
0166
0167       const double v_rot_y = fullAlignments[senId].getRotY();
0168       const double u_rot_y = (fullAlignments[senId].getRotYUnc() > 0.) ? fullAlignments[senId].getRotYUnc() : 1.;
0169       const double w_rot_y = (equalWeights) ? 1. : 1. / u_rot_y;
0170       sw2_rot_y += w_rot_y * w_rot_y;
0171       svw2_rot_y += v_rot_y * w_rot_y * w_rot_y;
0172       su2w4_rot_y += u_rot_y * u_rot_y * w_rot_y * w_rot_y * w_rot_y * w_rot_y;
0173
0174       const double v_rot_z = fullAlignments[senId].getRotZ();
0175       const double u_rot_z = (fullAlignments[senId].getRotZUnc() > 0.) ? fullAlignments[senId].getRotZUnc() : 1.;
0176       const double w_rot_z = (equalWeights) ? 1. : 1. / u_rot_z;
0177       sw2_rot_z += w_rot_z * w_rot_z;
0178       svw2_rot_z += v_rot_z * w_rot_z * w_rot_z;
0179       su2w4_rot_z += u_rot_z * u_rot_z * w_rot_z * w_rot_z * w_rot_z * w_rot_z;
0180
0181       if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0182         auto d2 = geometry.get(senId).getDirectionData(2);
0183
0184         B(idx, 0) = d2.dx;
0185         B(idx, 1) = d2.dy;
0186
0187         M(idx) = d2.dx * fullAlignments[senId].getShX() + d2.dy * fullAlignments[senId].getShY();
0188         double unc =
0189             sqrt(pow(d2.dx * fullAlignments[senId].getShXUnc(), 2) + pow(d2.dy * fullAlignments[senId].getShYUnc(), 2));
0190         if (unc <= 0.)
0191           unc = 1.;
0192
0193         Vi(idx, idx) = (equalWeights) ? 1. : 1. / unc / unc;
0194         VarM(idx, idx) = unc * unc;
0195
0196         idx += 1;
0197       }
0198
0199       if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0200         B(idx + 0, 0) = 1.;
0201         B(idx + 0, 1) = 0.;
0202         M(idx + 0) = fullAlignments[senId].getShX();
0203         double x_unc = fullAlignments[senId].getShXUnc();
0204         if (x_unc <= 0.)
0205           x_unc = 1.;
0206         Vi(idx + 0, idx + 0) = (equalWeights) ? 1. : 1. / x_unc / x_unc;
0207         VarM(idx + 0, idx + 0) = x_unc * x_unc;
0208
0209         B(idx + 1, 0) = 0.;
0210         B(idx + 1, 1) = 1.;
0211         M(idx + 1) = fullAlignments[senId].getShY();
0212         double y_unc = fullAlignments[senId].getShYUnc();
0213         if (y_unc <= 0.)
0214           y_unc = 1.;
0215         Vi(idx + 1, idx + 1) = (equalWeights) ? 1. : 1. / y_unc / y_unc;
0216         VarM(idx + 1, idx + 1) = y_unc * y_unc;
0217
0218         idx += 2;
0219       }
0220     }
0221
0222     // calculate per-RP alignment
0223     TMatrixD BT(TMatrixD::kTransposed, B);
0224     TMatrixD BTViB(BT, TMatrixD::kMult, Vi * B);
0225     TMatrixD BTViBi(TMatrixD::kInverted, BTViB);
0226     TMatrixD S(BTViBi * BT * Vi);
0227     TMatrixD ST(TMatrixD::kTransposed, S);
0228     TVectorD th_B(2);
0229     th_B = S * M;
0230     TMatrixD VarTh_B(S * VarM * ST);
0231
0232     const double m_sh_x = th_B[0], m_sh_x_unc = sqrt(VarTh_B(0, 0));
0233     const double m_sh_y = th_B[1], m_sh_y_unc = sqrt(VarTh_B(1, 1));
0234     const double m_sh_z = svw2_sh_z / sw2_sh_z, m_sh_z_unc = sqrt(su2w4_sh_z) / sw2_sh_z;
0235
0236     const double m_rot_x = svw2_rot_x / sw2_rot_x, m_rot_x_unc = sqrt(su2w4_rot_x) / sw2_rot_x;
0237     const double m_rot_y = svw2_rot_y / sw2_rot_y, m_rot_y_unc = sqrt(su2w4_rot_y) / sw2_rot_y;
0238     const double m_rot_z = svw2_rot_z / sw2_rot_z, m_rot_z_unc = sqrt(su2w4_rot_z) / sw2_rot_z;
0239
0240     if (verbosity) {
0241       printf("    m_sh_x = (%.1f +- %.1f) um, m_sh_y = (%.1f +- %.1f) um, m_sh_z = (%.1f +- %.1f) mm\n",
0242              m_sh_x * 1E3,
0243              m_sh_x_unc * 1E3,
0244              m_sh_y * 1E3,
0245              m_sh_y_unc * 1E3,
0246              m_sh_z,
0247              m_sh_z_unc);
0248       printf("    m_rot_x = (%.1f +- %.1f) mrad, m_rot_y = (%.1f +- %.1f)  mrad, m_rot_z = (%.1f +- %.1f) mrad\n",
0249              m_rot_x * 1E3,
0250              m_rot_x_unc * 1E3,
0251              m_rot_y * 1E3,
0252              m_rot_y_unc * 1E3,
0253              m_rot_z * 1E3,
0254              m_rot_z_unc * 1E3);
0255     }
0256
0257     factoredAlignments.addRPCorrection(rpId,
0258                                        CTPPSRPAlignmentCorrectionData(m_sh_x,
0259                                                                       m_sh_x_unc,
0260                                                                       m_sh_y,
0261                                                                       m_sh_y_unc,
0262                                                                       m_sh_z,
0263                                                                       m_sh_z_unc,
0264                                                                       m_rot_x,
0265                                                                       m_rot_x_unc,
0266                                                                       m_rot_y,
0267                                                                       m_rot_y_unc,
0268                                                                       m_rot_z,
0269                                                                       m_rot_z_unc));
0270
0271     // calculate residuals
0272     for (const auto &senId : sensors) {
0273       CTPPSRPAlignmentCorrectionData rc;
0274
0275       if (rpId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0276         auto d2 = geometry.get(senId).getDirectionData(2);
0277
0278         const double de_s =
0279             d2.dx * (fullAlignments[senId].getShX() - m_sh_x) + d2.dy * (fullAlignments[senId].getShY() - m_sh_y);
0280         const double de_s_unc =
0281             std::abs(d2.dx * fullAlignments[senId].getShXUnc() +
0282                      d2.dy * fullAlignments[senId].getShYUnc());  // the x and y components are fully correlated
0283
0284         rc = CTPPSRPAlignmentCorrectionData(d2.dx * de_s,
0285                                             d2.dx * de_s_unc,
0286                                             d2.dy * de_s,
0287                                             d2.dy * de_s_unc,
0288                                             fullAlignments[senId].getShZ() - m_sh_z,
0289                                             fullAlignments[senId].getShZUnc(),
0290                                             fullAlignments[senId].getRotX() - m_rot_x,
0291                                             fullAlignments[senId].getRotXUnc(),
0292                                             fullAlignments[senId].getRotY() - m_rot_y,
0293                                             fullAlignments[senId].getRotYUnc(),
0294                                             fullAlignments[senId].getRotZ() - m_rot_z,
0295                                             fullAlignments[senId].getRotZUnc());
0296       }
0297
0298       if (rpId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0299         rc = CTPPSRPAlignmentCorrectionData(fullAlignments[senId].getShX() - m_sh_x,
0300                                             fullAlignments[senId].getShXUnc(),
0301                                             fullAlignments[senId].getShY() - m_sh_y,
0302                                             fullAlignments[senId].getShYUnc(),
0303                                             fullAlignments[senId].getShZ() - m_sh_z,
0304                                             fullAlignments[senId].getShZUnc(),
0305                                             fullAlignments[senId].getRotX() - m_rot_x,
0306                                             fullAlignments[senId].getRotXUnc(),
0307                                             fullAlignments[senId].getRotY() - m_rot_y,
0308                                             fullAlignments[senId].getRotYUnc(),
0309                                             fullAlignments[senId].getRotZ() - m_rot_z,
0310                                             fullAlignments[senId].getRotZUnc());
0311       }
0312
0313       factoredAlignments.addSensorCorrection(senId, rc);
0314     }
0315   }
0316 }