File indexing completed on 2024-04-06 11:56:43
0001 #include <iostream>
0002
0003 #include "CLHEP/Matrix/Matrix.h"
0004 #include "TMatrixDEigen.h"
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "Alignment/MuonAlignmentAlgorithms/plugins/CSCChamberFitter.h"
0009
0010 const double infinity =
0011 0.1;
0012
0013
0014 CSCChamberFitter::CSCChamberFitter(const edm::ParameterSet &iConfig,
0015 std::vector<CSCPairResidualsConstraint *> &residualsConstraints) {
0016 m_name = iConfig.getParameter<std::string>("name");
0017 m_alignables = iConfig.getParameter<std::vector<std::string> >("alignables");
0018 if (m_alignables.empty()) {
0019 throw cms::Exception("BadConfig") << "Fitter " << m_name << " has no alignables!" << std::endl;
0020 }
0021
0022 int i = 0;
0023 for (std::vector<std::string>::const_iterator alignable = m_alignables.begin(); alignable != m_alignables.end();
0024 ++alignable) {
0025 if (alignableId(*alignable) == -1)
0026 m_frames.push_back(i);
0027 i++;
0028 }
0029
0030 m_fixed = -1;
0031 std::string fixed = iConfig.getParameter<std::string>("fixed");
0032 if (!fixed.empty()) {
0033 int i = 0;
0034 for (std::vector<std::string>::const_iterator alignable = m_alignables.begin(); alignable != m_alignables.end();
0035 ++alignable) {
0036 if (fixed == *alignable) {
0037 m_fixed = i;
0038 }
0039 i++;
0040 }
0041 if (m_fixed == -1)
0042 throw cms::Exception("BadConfig") << "Cannot fix unrecognized alignable " << fixed << std::endl;
0043 }
0044
0045 int numConstraints = 0;
0046 std::vector<edm::ParameterSet> constraints = iConfig.getParameter<std::vector<edm::ParameterSet> >("constraints");
0047 for (std::vector<edm::ParameterSet>::const_iterator constraint = constraints.begin(); constraint != constraints.end();
0048 ++constraint) {
0049 int i = index(constraint->getParameter<std::string>("i"));
0050 int j = index(constraint->getParameter<std::string>("j"));
0051 double value = constraint->getParameter<double>("value");
0052 double error = constraint->getParameter<double>("error");
0053
0054 if (i < 0)
0055 throw cms::Exception("BadConfig") << "Unrecognized alignable " << constraint->getParameter<std::string>("i")
0056 << " in constraint " << numConstraints << " of fitter " << m_name << std::endl;
0057 if (j < 0)
0058 throw cms::Exception("BadConfig") << "Unrecognized alignable " << constraint->getParameter<std::string>("j")
0059 << " in constraint " << numConstraints << " of fitter " << m_name << std::endl;
0060 if (error <= 0.)
0061 throw cms::Exception("BadConfig") << "Non-positive uncertainty in constraint " << numConstraints << " of fitter "
0062 << m_name << std::endl;
0063 if (i == j)
0064 throw cms::Exception("BadConfig") << "Self-connection from " << constraint->getParameter<std::string>("i")
0065 << " to " << constraint->getParameter<std::string>("j")
0066 << " is not allowed in constraint " << numConstraints << " of fitter " << m_name
0067 << std::endl;
0068
0069 m_constraints.push_back(new CSCPairConstraint(i, j, value, error));
0070 numConstraints++;
0071 }
0072
0073
0074 for (unsigned int i = 0; i < m_alignables.size(); i++) {
0075 std::string alignable_i = m_alignables[i];
0076 long id_i = alignableId(alignable_i);
0077 if (id_i != -1) {
0078 CSCDetId cscid_i(id_i);
0079
0080 for (unsigned int j = 0; j < m_alignables.size(); j++) {
0081 std::string alignable_j = m_alignables[j];
0082 long id_j = alignableId(alignable_j);
0083 if (i != j && id_j != -1) {
0084 CSCDetId cscid_j(id_j);
0085
0086 if (!(cscid_i.station() == 1 && cscid_i.ring() == 3 && cscid_j.station() == 1 && cscid_j.ring() == 3)) {
0087 int next_chamber = cscid_i.chamber() + 1;
0088 if (cscid_i.station() > 1 && cscid_i.ring() == 1 && next_chamber == 19)
0089 next_chamber = 1;
0090 else if (!(cscid_i.station() > 1 && cscid_i.ring() == 1) && next_chamber == 37)
0091 next_chamber = 1;
0092 if (cscid_i.endcap() == cscid_j.endcap() && cscid_i.station() == cscid_j.station() &&
0093 cscid_i.ring() == cscid_j.ring() && next_chamber == cscid_j.chamber()) {
0094 CSCPairResidualsConstraint *residualsConstraint =
0095 new CSCPairResidualsConstraint(residualsConstraints.size(), i, j, cscid_i, cscid_j);
0096 m_constraints.push_back(residualsConstraint);
0097 residualsConstraints.push_back(residualsConstraint);
0098 numConstraints++;
0099 }
0100 }
0101 }
0102 }
0103 }
0104 }
0105
0106 std::map<int, bool> touched;
0107 for (unsigned int i = 0; i < m_alignables.size(); i++)
0108 touched[i] = false;
0109 walk(touched, 0);
0110 for (unsigned int i = 0; i < m_alignables.size(); i++) {
0111 if (!touched[i])
0112 throw cms::Exception("BadConfig") << "Fitter " << m_name << " is not a connected graph (no way to get to "
0113 << m_alignables[i] << " from " << m_alignables[0] << ", for instance)"
0114 << std::endl;
0115 }
0116 }
0117
0118 int CSCChamberFitter::index(std::string alignable) const {
0119 int i = 0;
0120 for (std::vector<std::string>::const_iterator a = m_alignables.begin(); a != m_alignables.end(); ++a) {
0121 if (*a == alignable)
0122 return i;
0123 i++;
0124 }
0125 return -1;
0126 }
0127
0128 void CSCChamberFitter::walk(std::map<int, bool> &touched, int alignable) const {
0129 touched[alignable] = true;
0130
0131 for (std::vector<CSCPairConstraint *>::const_iterator constraint = m_constraints.begin();
0132 constraint != m_constraints.end();
0133 ++constraint) {
0134 if (alignable == (*constraint)->i() || alignable == (*constraint)->j()) {
0135 if (!touched[(*constraint)->i()])
0136 walk(touched, (*constraint)->i());
0137 if (!touched[(*constraint)->j()])
0138 walk(touched, (*constraint)->j());
0139 }
0140 }
0141 }
0142
0143 long CSCChamberFitter::alignableId(std::string alignable) const {
0144 if (alignable.size() != 9)
0145 return -1;
0146
0147 if (alignable[0] == 'M' && alignable[1] == 'E') {
0148 int endcap = -1;
0149 if (alignable[2] == '+')
0150 endcap = 1;
0151 else if (alignable[2] == '-')
0152 endcap = 2;
0153
0154 if (endcap != -1) {
0155 int station = -1;
0156 if (alignable[3] == '1')
0157 station = 1;
0158 else if (alignable[3] == '2')
0159 station = 2;
0160 else if (alignable[3] == '3')
0161 station = 3;
0162 else if (alignable[3] == '4')
0163 station = 4;
0164
0165 if (alignable[4] == '/' && station != -1) {
0166 int ring = -1;
0167 if (alignable[5] == '1')
0168 ring = 1;
0169 else if (alignable[5] == '2')
0170 ring = 2;
0171 else if (alignable[5] == '3')
0172 ring = 3;
0173 else if (alignable[5] == '4')
0174 ring = 4;
0175 if (station > 1 && ring > 2)
0176 return -1;
0177
0178 if (alignable[6] == '/' && ring != -1) {
0179 int chamber = -1;
0180 if (alignable[7] == '0' && alignable[8] == '1')
0181 chamber = 1;
0182 else if (alignable[7] == '0' && alignable[8] == '2')
0183 chamber = 2;
0184 else if (alignable[7] == '0' && alignable[8] == '3')
0185 chamber = 3;
0186 else if (alignable[7] == '0' && alignable[8] == '4')
0187 chamber = 4;
0188 else if (alignable[7] == '0' && alignable[8] == '5')
0189 chamber = 5;
0190 else if (alignable[7] == '0' && alignable[8] == '6')
0191 chamber = 6;
0192 else if (alignable[7] == '0' && alignable[8] == '7')
0193 chamber = 7;
0194 else if (alignable[7] == '0' && alignable[8] == '8')
0195 chamber = 8;
0196 else if (alignable[7] == '0' && alignable[8] == '9')
0197 chamber = 9;
0198 else if (alignable[7] == '1' && alignable[8] == '0')
0199 chamber = 10;
0200 else if (alignable[7] == '1' && alignable[8] == '1')
0201 chamber = 11;
0202 else if (alignable[7] == '1' && alignable[8] == '2')
0203 chamber = 12;
0204 else if (alignable[7] == '1' && alignable[8] == '3')
0205 chamber = 13;
0206 else if (alignable[7] == '1' && alignable[8] == '4')
0207 chamber = 14;
0208 else if (alignable[7] == '1' && alignable[8] == '5')
0209 chamber = 15;
0210 else if (alignable[7] == '1' && alignable[8] == '6')
0211 chamber = 16;
0212 else if (alignable[7] == '1' && alignable[8] == '7')
0213 chamber = 17;
0214 else if (alignable[7] == '1' && alignable[8] == '8')
0215 chamber = 18;
0216 else if (alignable[7] == '1' && alignable[8] == '9')
0217 chamber = 19;
0218 else if (alignable[7] == '2' && alignable[8] == '0')
0219 chamber = 20;
0220 else if (alignable[7] == '2' && alignable[8] == '1')
0221 chamber = 21;
0222 else if (alignable[7] == '2' && alignable[8] == '2')
0223 chamber = 22;
0224 else if (alignable[7] == '2' && alignable[8] == '3')
0225 chamber = 23;
0226 else if (alignable[7] == '2' && alignable[8] == '4')
0227 chamber = 24;
0228 else if (alignable[7] == '2' && alignable[8] == '5')
0229 chamber = 25;
0230 else if (alignable[7] == '2' && alignable[8] == '6')
0231 chamber = 26;
0232 else if (alignable[7] == '2' && alignable[8] == '7')
0233 chamber = 27;
0234 else if (alignable[7] == '2' && alignable[8] == '8')
0235 chamber = 28;
0236 else if (alignable[7] == '2' && alignable[8] == '9')
0237 chamber = 29;
0238 else if (alignable[7] == '3' && alignable[8] == '0')
0239 chamber = 30;
0240 else if (alignable[7] == '3' && alignable[8] == '1')
0241 chamber = 31;
0242 else if (alignable[7] == '3' && alignable[8] == '2')
0243 chamber = 32;
0244 else if (alignable[7] == '3' && alignable[8] == '3')
0245 chamber = 33;
0246 else if (alignable[7] == '3' && alignable[8] == '4')
0247 chamber = 34;
0248 else if (alignable[7] == '3' && alignable[8] == '5')
0249 chamber = 35;
0250 else if (alignable[7] == '3' && alignable[8] == '6')
0251 chamber = 36;
0252
0253 if (station > 1 && ring == 1 && chamber > 18)
0254 return -1;
0255
0256 if (chamber != -1) {
0257 return CSCDetId(endcap, station, ring, chamber, 0).rawId();
0258 }
0259 }
0260 }
0261 }
0262 }
0263
0264 return -1;
0265 }
0266
0267 bool CSCChamberFitter::isFrame(int i) const {
0268 for (std::vector<int>::const_iterator frame = m_frames.begin(); frame != m_frames.end(); ++frame) {
0269 if (i == *frame)
0270 return true;
0271 }
0272 return false;
0273 }
0274
0275 double CSCChamberFitter::chi2(const AlgebraicVector &A, double lambda) const {
0276 double sumFixed = 0.;
0277
0278 if (m_fixed == -1) {
0279 for (unsigned int i = 0; i < m_alignables.size(); i++) {
0280 if (!isFrame(i)) {
0281 sumFixed += A[i];
0282 }
0283 }
0284 } else {
0285 sumFixed = A[m_fixed];
0286 }
0287
0288 double s = lambda * sumFixed * sumFixed;
0289 for (std::vector<CSCPairConstraint *>::const_iterator constraint = m_constraints.begin();
0290 constraint != m_constraints.end();
0291 ++constraint) {
0292 if ((*constraint)->valid()) {
0293 s += pow((*constraint)->value() - A[(*constraint)->i()] + A[(*constraint)->j()], 2) / (*constraint)->error() /
0294 (*constraint)->error();
0295 }
0296 }
0297 return s;
0298 }
0299
0300 double CSCChamberFitter::lhsVector(int k) const {
0301 double s = 0.;
0302 for (std::vector<CSCPairConstraint *>::const_iterator constraint = m_constraints.begin();
0303 constraint != m_constraints.end();
0304 ++constraint) {
0305 if ((*constraint)->valid()) {
0306 double d = 2. * (*constraint)->value() / (*constraint)->error() / (*constraint)->error();
0307 if ((*constraint)->i() == k)
0308 s += d;
0309 if ((*constraint)->j() == k)
0310 s -= d;
0311 }
0312 }
0313 return s;
0314 }
0315
0316 double CSCChamberFitter::hessian(int k, int l, double lambda) const {
0317 double s = 0.;
0318
0319 if (m_fixed == -1) {
0320 if (!isFrame(k) && !isFrame(l))
0321 s += 2. * lambda;
0322 } else {
0323 if (k == l && l == m_fixed)
0324 s += 2. * lambda;
0325 }
0326
0327 for (std::vector<CSCPairConstraint *>::const_iterator constraint = m_constraints.begin();
0328 constraint != m_constraints.end();
0329 ++constraint) {
0330 double d = 2. / infinity / infinity;
0331 if ((*constraint)->valid()) {
0332 d = 2. / (*constraint)->error() / (*constraint)->error();
0333 }
0334
0335 if (k == l && ((*constraint)->i() == k || (*constraint)->j() == k))
0336 s += d;
0337 if (((*constraint)->i() == k && (*constraint)->j() == l) || ((*constraint)->j() == k && (*constraint)->i() == l))
0338 s -= d;
0339 }
0340 return s;
0341 }
0342
0343 bool CSCChamberFitter::fit(std::vector<CSCAlignmentCorrections *> &corrections) const {
0344 double lambda = 1. / infinity / infinity;
0345
0346 AlgebraicVector A(m_alignables.size());
0347 AlgebraicVector V(m_alignables.size());
0348 AlgebraicMatrix M(m_alignables.size(), m_alignables.size());
0349
0350 for (unsigned int k = 0; k < m_alignables.size(); k++) {
0351 A[k] = 0.;
0352 V[k] = lhsVector(k);
0353 for (unsigned int l = 0; l < m_alignables.size(); l++) {
0354 M[k][l] = hessian(k, l, lambda);
0355 }
0356 }
0357
0358 double oldchi2 = chi2(A, lambda);
0359
0360 int ierr;
0361 M.invert(ierr);
0362 if (ierr != 0) {
0363 edm::LogError("CSCOverlapsAlignmentAlgorithm")
0364 << "Matrix inversion failed for fitter " << m_name << " matrix is " << M << std::endl;
0365 return false;
0366 }
0367
0368 A = M * V;
0369
0370
0371 CSCAlignmentCorrections *correction = new CSCAlignmentCorrections(m_name, oldchi2, chi2(A, lambda));
0372
0373 for (unsigned int i = 0; i < m_alignables.size(); i++) {
0374 if (!isFrame(i)) {
0375 correction->insertCorrection(m_alignables[i], CSCDetId(alignableId(m_alignables[i])), A[i]);
0376 }
0377 }
0378
0379
0380
0381
0382 TMatrixD tmatrix(m_alignables.size(), m_alignables.size());
0383 for (unsigned int i = 0; i < m_alignables.size(); i++) {
0384 for (unsigned int j = 0; j < m_alignables.size(); j++) {
0385 tmatrix[i][j] = M[i][j];
0386 }
0387 }
0388 TMatrixDEigen tmatrixdeigen(tmatrix);
0389 const TMatrixD &basis = tmatrixdeigen.GetEigenVectors();
0390 TMatrixD invbasis = tmatrixdeigen.GetEigenVectors();
0391 invbasis.Invert();
0392 TMatrixD diagonalized = invbasis * (tmatrix * basis);
0393
0394 for (unsigned int i = 0; i < m_alignables.size(); i++) {
0395 std::vector<double> coefficient;
0396 std::vector<std::string> modename;
0397 std::vector<long> modeid;
0398 for (unsigned int j = 0; j < m_alignables.size(); j++) {
0399 coefficient.push_back(invbasis[i][j]);
0400 modename.push_back(m_alignables[j]);
0401 modeid.push_back(alignableId(m_alignables[j]));
0402 }
0403
0404 correction->insertMode(
0405 coefficient, modename, modeid, sqrt(2. * fabs(diagonalized[i][i])) * (diagonalized[i][i] >= 0. ? 1. : -1.));
0406 }
0407
0408 for (std::vector<CSCPairConstraint *>::const_iterator constraint = m_constraints.begin();
0409 constraint != m_constraints.end();
0410 ++constraint) {
0411 if ((*constraint)->valid()) {
0412 double residual = (*constraint)->value() - A[(*constraint)->i()] + A[(*constraint)->j()];
0413 correction->insertResidual(m_alignables[(*constraint)->i()],
0414 m_alignables[(*constraint)->j()],
0415 (*constraint)->value(),
0416 (*constraint)->error(),
0417 residual,
0418 residual / (*constraint)->error());
0419 }
0420 }
0421
0422 corrections.push_back(correction);
0423 return true;
0424 }
0425
0426 void CSCChamberFitter::radiusCorrection(AlignableNavigator *alignableNavigator,
0427 AlignmentParameterStore *alignmentParameterStore,
0428 bool combineME11) const {
0429 double sum_phipos_residuals = 0.;
0430 double num_valid = 0.;
0431 double sum_radius = 0.;
0432 double num_total = 0.;
0433 for (std::vector<CSCPairConstraint *>::const_iterator constraint = m_constraints.begin();
0434 constraint != m_constraints.end();
0435 ++constraint) {
0436 CSCPairResidualsConstraint *residualsConstraint = dynamic_cast<CSCPairResidualsConstraint *>(*constraint);
0437 if (residualsConstraint != nullptr) {
0438 if (residualsConstraint->valid()) {
0439 sum_phipos_residuals += residualsConstraint->value();
0440 num_valid += 1.;
0441 }
0442
0443 sum_radius += residualsConstraint->radius(true);
0444 num_total += 1.;
0445 }
0446 }
0447 if (num_valid == 0. || num_total == 0.)
0448 return;
0449 double average_phi_residual = sum_phipos_residuals / num_valid;
0450 double average_radius = sum_radius / num_total;
0451
0452 double radial_correction = average_phi_residual * average_radius * num_total / (2. * M_PI);
0453
0454 for (std::vector<CSCPairConstraint *>::const_iterator constraint = m_constraints.begin();
0455 constraint != m_constraints.end();
0456 ++constraint) {
0457 CSCPairResidualsConstraint *residualsConstraint = dynamic_cast<CSCPairResidualsConstraint *>(*constraint);
0458 if (residualsConstraint != nullptr) {
0459 const DetId id(residualsConstraint->id_i());
0460 Alignable *alignable = alignableNavigator->alignableFromDetId(id).alignable();
0461 Alignable *also = nullptr;
0462 if (combineME11 && residualsConstraint->id_i().station() == 1 && residualsConstraint->id_i().ring() == 1) {
0463 CSCDetId alsoid(residualsConstraint->id_i().endcap(), 1, 4, residualsConstraint->id_i().chamber(), 0);
0464 const DetId alsoid2(alsoid);
0465 also = alignableNavigator->alignableFromDetId(alsoid2).alignable();
0466 }
0467
0468 AlgebraicVector params(6);
0469 AlgebraicSymMatrix cov(6);
0470
0471 params[1] = radial_correction;
0472 cov[1][1] = 1e-6;
0473
0474 AlignmentParameters *parnew = alignable->alignmentParameters()->cloneFromSelected(params, cov);
0475 alignable->setAlignmentParameters(parnew);
0476 alignmentParameterStore->applyParameters(alignable);
0477 alignable->alignmentParameters()->setValid(true);
0478 if (also != nullptr) {
0479 AlignmentParameters *parnew2 = also->alignmentParameters()->cloneFromSelected(params, cov);
0480 also->setAlignmentParameters(parnew2);
0481 alignmentParameterStore->applyParameters(also);
0482 also->alignmentParameters()->setValid(true);
0483 }
0484 }
0485 }
0486 }