AlignmentRelative/src/LocalTrackFitter.cc

0001 /****************************************************************************
0002 * Authors:
0003 *  Jan Kašpar (jan.kaspar@gmail.com)
0004 ****************************************************************************/
0005
0006 #include "CalibPPS/AlignmentRelative/interface/LocalTrackFitter.h"
0007
0008 #include "CalibPPS/AlignmentRelative/interface/Utilities.h"
0009
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011
0012 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0013 #include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h"
0014
0015 #include <set>
0016
0017 #include "TMatrixD.h"
0018 #include "TVectorD.h"
0019
0020 using namespace std;
0021
0022 //----------------------------------------------------------------------------------------------------
0023
0024 LocalTrackFitter::LocalTrackFitter(const edm::ParameterSet &ps)
0025     : verbosity(ps.getUntrackedParameter<unsigned int>("verbosity", 0)),
0026       minimumHitsPerProjectionPerRP(ps.getParameter<unsigned int>("minimumHitsPerProjectionPerRP")),
0027       maxResidualToSigma(ps.getParameter<double>("maxResidualToSigma")) {}
0028
0029 //----------------------------------------------------------------------------------------------------
0030
0031 bool LocalTrackFitter::fit(HitCollection &selection, const AlignmentGeometry &geometry, LocalTrackFit &trackFit) const {
0032   if (verbosity > 5)
0033     printf(">> LocalTrackFitter::Fit\n");
0034
0035   bool selectionChanged = true;
0036   unsigned int loopCounter = 0;
0037   while (selectionChanged) {
0038     // fit/outlier-removal loop
0039     while (selectionChanged) {
0040       if (verbosity > 5)
0041         printf("* fit loop %u\n", loopCounter++);
0042
0043       bool fitFailed = false;
0044       fitAndRemoveOutliers(selection, geometry, trackFit, fitFailed, selectionChanged);
0045
0046       if (fitFailed) {
0047         if (verbosity > 5)
0048           printf("\tFIT FAILED\n");
0049         return false;
0050       }
0051     }
0052
0053     // remove pots with too few active planes
0054     removeInsufficientPots(selection, selectionChanged);
0055   }
0056
0057   return true;
0058 }
0059
0060 //----------------------------------------------------------------------------------------------------
0061
0062 void LocalTrackFitter::fitAndRemoveOutliers(HitCollection &selection,
0063                                             const AlignmentGeometry &geometry,
0064                                             LocalTrackFit &trackFit,
0065                                             bool &failed,
0066                                             bool &selectionChanged) const {
0067   if (verbosity > 5)
0068     printf(" - LocalTrackFitter::FitAndRemoveOutliers\n");
0069
0070   if (selection.empty()) {
0071     failed = true;
0072     return;
0073   }
0074
0075   // check if input size is sufficient
0076   if (selection.size() < 4) {
0077     failed = true;
0078     return;
0079   }
0080
0081   // build matrices and vectors
0082   TMatrixD A(selection.size(), 4);
0083   TMatrixD Vi(selection.size(), selection.size());
0084   TVectorD measVec(selection.size());
0085   unsigned int j = 0;
0086   for (auto it = selection.begin(); it != selection.end(); ++it, ++j) {
0087     const unsigned int &detId = it->id;
0088
0089     const DetGeometry &d = geometry.get(detId);
0090     const auto &dirData = d.getDirectionData(it->dirIdx);
0091
0092     A(j, 0) = it->z * dirData.dx;
0093     A(j, 1) = dirData.dx;
0094     A(j, 2) = it->z * dirData.dy;
0095     A(j, 3) = dirData.dy;
0096
0097     measVec(j) = it->position + dirData.s - (it->z - d.z) * dirData.dz;  // in mm
0098
0099     Vi(j, j) = 1. / it->sigma / it->sigma;
0100   }
0101
0102   // evaluate local track parameter estimates (h stands for hat)
0103   TMatrixD AT(4, selection.size());
0104   AT.Transpose(A);
0105   TMatrixD ATViA(4, 4);
0106   ATViA = AT * Vi * A;
0107   TMatrixD ATViAI(ATViA);
0108
0109   try {
0110     ATViAI = ATViA.Invert();
0111   } catch (cms::Exception &e) {
0112     failed = true;
0113     return;
0114   }
0115
0116   TVectorD theta(4);
0117   theta = ATViAI * AT * Vi * measVec;
0118
0119   // residuals
0120   TVectorD R(measVec);
0121   R -= A * theta;
0122
0123   // save results to trackFit
0124   trackFit.ax = theta(0);
0125   trackFit.bx = theta(1);
0126   trackFit.ay = theta(2);
0127   trackFit.by = theta(3);
0128   trackFit.z0 = geometry.z0;
0129   trackFit.ndf = selection.size() - 4;
0130   trackFit.chi_sq = 0;
0131   for (int i = 0; i < R.GetNrows(); i++)
0132     trackFit.chi_sq += R(i) * R(i) * Vi(i, i);
0133
0134   if (verbosity > 5) {
0135     printf("    ax = %.3f mrad, bx = %.4f mm, ay = %.3f mrad, by = %.4f mm, z0 = %.3f mm\n",
0136            trackFit.ax * 1E3,
0137            trackFit.bx,
0138            trackFit.ay * 1E3,
0139            trackFit.by,
0140            trackFit.z0);
0141     printf("    ndof = %i, chi^2/ndof/si^2 = %.3f\n", trackFit.ndf, trackFit.chi_sq / trackFit.ndf);
0142   }
0143
0144   // check residuals
0145   selectionChanged = false;
0146   TVectorD interpolation(A * theta);
0147   j = 0;
0148   for (auto it = selection.begin(); it != selection.end(); ++j) {
0149     if (verbosity > 5) {
0150       printf("        %2u, ", j);
0151       printId(it->id);
0152       printf(", dirIdx=%u: interpol = %+8.1f um, residual = %+6.1f um, residual / sigma = %+6.2f\n",
0153              it->dirIdx,
0154              interpolation[j] * 1E3,
0155              R[j] * 1E3,
0156              R[j] / it->sigma);
0157     }
0158
0159     double resToSigma = R[j] / it->sigma;
0160     if (fabs(resToSigma) > maxResidualToSigma) {
0161       selection.erase(it);
0162       selectionChanged = true;
0163       if (verbosity > 5)
0164         printf("            Removed\n");
0165     } else
0166       ++it;
0167   }
0168 }
0169
0170 //----------------------------------------------------------------------------------------------------
0171
0172 void LocalTrackFitter::removeInsufficientPots(HitCollection &selection, bool &selectionChanged) const {
0173   if (verbosity > 5)
0174     printf(" - RemoveInsufficientPots\n");
0175
0176   selectionChanged = false;
0177
0178   // map: RP id -> (active planes in projection 1, active planes in projection 2)
0179   map<unsigned int, pair<set<unsigned int>, set<unsigned int> > > planeMap;
0180   for (auto it = selection.begin(); it != selection.end(); ++it) {
0181     CTPPSDetId senId(it->id);
0182     const unsigned int rpId = senId.rpId();
0183
0184     if (senId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0185       const unsigned int plane = TotemRPDetId(it->id).plane();
0186       if ((plane % 2) == 0)
0187         planeMap[rpId].first.insert(senId);
0188       else
0189         planeMap[rpId].second.insert(senId);
0190     }
0191
0192     if (senId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0193       planeMap[rpId].first.insert(senId);
0194       planeMap[rpId].second.insert(senId);
0195     }
0196   }
0197
0198   // remove RPs with insufficient information
0199   selectionChanged = false;
0200
0201   for (auto it = planeMap.begin(); it != planeMap.end(); ++it) {
0202     if (it->second.first.size() < minimumHitsPerProjectionPerRP ||
0203         it->second.second.size() < minimumHitsPerProjectionPerRP) {
0204       if (verbosity > 5)
0205         printf("\tRP %u: projection1 = %lu, projection 2 = %lu\n",
0206                it->first,
0207                it->second.first.size(),
0208                it->second.second.size());
0209
0210       // remove all hits from that RP
0211       for (auto dit = selection.begin(); dit != selection.end();) {
0212         if (it->first == CTPPSDetId(dit->id).rpId()) {
0213           if (verbosity > 5)
0214             printf("\t\tremoving %u\n", dit->id);
0215           selection.erase(dit);
0216           selectionChanged = true;
0217         } else {
0218           dit++;
0219         }
0220       }
0221     }
0222   }
0223 }