File indexing completed on 2024-04-06 11:58:34
0001
0002
0003
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
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
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
0076 if (selection.size() < 4) {
0077 failed = true;
0078 return;
0079 }
0080
0081
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;
0098
0099 Vi(j, j) = 1. / it->sigma / it->sigma;
0100 }
0101
0102
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
0120 TVectorD R(measVec);
0121 R -= A * theta;
0122
0123
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
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
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
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
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 }