Back to home page

Project CMSSW displayed by LXR
 
 

    

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;  // this is huge because all alignments are angles in radians; but we need a not-too-large value for numerical stability
0012           // should become a parameter someday
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   // insert CSCPairResidualsConstraints
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;  // that's the alignment step
0369 
0370   ///// everything else is for reporting
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   // we have to switch to a completely different linear algebrea
0380   // package because CLHEP doesn't compute
0381   // eigenvectors/diagonalization (?!?)
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 }
[ 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!