File indexing completed on 2023-03-17 10:42:30
0001
0002
0003
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
0100 expandedAlignments.clear();
0101 factoredAlignments.clear();
0102
0103
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
0110 fullAlignments[sensorId] = inputAlignments.getFullSensorCorrection(sensorId, false);
0111
0112 sensorsPerRP[CTPPSDetId(sensorId).rpId()].insert(sensorId);
0113 }
0114
0115
0116 for (const auto &it : fullAlignments) {
0117 expandedAlignments.setSensorCorrection(it.first, it.second);
0118 }
0119
0120
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
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
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
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
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());
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 }