File indexing completed on 2024-04-06 12:26:14
0001
0002
0003
0004
0005
0006
0007 #include "RecoLocalMuon/GEMSegment/plugins/MuonSegFit.h"
0008 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0009
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 bool MuonSegFit::fit(void) {
0014 if (fitdone())
0015 return fitdone_;
0016 short n = nhits();
0017 if (n < 2) {
0018 edm::LogVerbatim("MuonSegFit") << "[MuonSegFit::fit] - cannot fit just 1 hit!!";
0019 } else if (n == 2) {
0020 fit2();
0021 } else if (2 * n <= MaxHits2) {
0022 fitlsq();
0023 } else {
0024 edm::LogVerbatim("MuonSegFit") << "[MuonSegFit::fit] - cannot fit more than " << MaxHits2 / 2 << " hits!!";
0025 }
0026 return fitdone_;
0027 }
0028
0029 void MuonSegFit::fit2(void) {
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044 MuonRecHitContainer::const_iterator ih = hits_.begin();
0045
0046
0047
0048 LocalPoint h1pos = (*ih)->localPosition();
0049 ++ih;
0050 LocalPoint h2pos = (*ih)->localPosition();
0051
0052
0053
0054 float dz = h2pos.z() - h1pos.z();
0055 if (dz != 0.0) {
0056 uslope_ = (h2pos.x() - h1pos.x()) / dz;
0057 vslope_ = (h2pos.y() - h1pos.y()) / dz;
0058 }
0059
0060 float uintercept = (h1pos.x() * h2pos.z() - h2pos.x() * h1pos.z()) / dz;
0061 float vintercept = (h1pos.y() * h2pos.z() - h2pos.y() * h1pos.z()) / dz;
0062
0063
0064 intercept_ = LocalPoint(uintercept, vintercept, 0.);
0065
0066
0067 setOutFromIP();
0068
0069
0070 chi2_ = 0.;
0071 ndof_ = 0;
0072
0073 fitdone_ = true;
0074 }
0075
0076 void MuonSegFit::fitlsq(void) {
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138 SMatrix4 M;
0139 SVector4 B;
0140
0141 MuonRecHitContainer::const_iterator ih = hits_.begin();
0142
0143
0144 for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0145 LocalPoint lp = (*ih)->localPosition();
0146
0147 double u = lp.x();
0148 double v = lp.y();
0149 double z = lp.z();
0150
0151
0152 SMatrixSym2 IC;
0153
0154 IC(0, 0) = (*ih)->localPositionError().xx();
0155 IC(1, 1) = (*ih)->localPositionError().yy();
0156
0157
0158 IC(1, 0) = (*ih)->localPositionError().xy();
0159
0160
0161 #ifdef EDM_ML_DEBUG
0162 edm::LogVerbatim("MuonSegFitMatrixDetails")
0163 << "[MuonSegFit::fit] 2x2 covariance matrix for this GEMRecHit :: [[" << IC(0, 0) << ", " << IC(0, 1) << "]["
0164 << IC(1, 0) << "," << IC(1, 1) << "]]";
0165 #endif
0166
0167
0168 bool ok = IC.Invert();
0169 if (!ok) {
0170 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::fit] Failed to invert covariance matrix: \n" << IC;
0171
0172 }
0173
0174 M(0, 0) += IC(0, 0);
0175 M(0, 1) += IC(0, 1);
0176 M(0, 2) += IC(0, 0) * z;
0177 M(0, 3) += IC(0, 1) * z;
0178 B(0) += u * IC(0, 0) + v * IC(0, 1);
0179
0180 M(1, 0) += IC(1, 0);
0181 M(1, 1) += IC(1, 1);
0182 M(1, 2) += IC(1, 0) * z;
0183 M(1, 3) += IC(1, 1) * z;
0184 B(1) += u * IC(1, 0) + v * IC(1, 1);
0185
0186 M(2, 0) += IC(0, 0) * z;
0187 M(2, 1) += IC(0, 1) * z;
0188 M(2, 2) += IC(0, 0) * z * z;
0189 M(2, 3) += IC(0, 1) * z * z;
0190 B(2) += (u * IC(0, 0) + v * IC(0, 1)) * z;
0191
0192 M(3, 0) += IC(1, 0) * z;
0193 M(3, 1) += IC(1, 1) * z;
0194 M(3, 2) += IC(1, 0) * z * z;
0195 M(3, 3) += IC(1, 1) * z * z;
0196 B(3) += (u * IC(1, 0) + v * IC(1, 1)) * z;
0197 }
0198
0199 SVector4 p;
0200 bool ok = M.Invert();
0201 if (!ok) {
0202 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::fit] Failed to invert matrix: \n" << M;
0203
0204 } else {
0205 p = M * B;
0206 }
0207
0208 #ifdef EDM_ML_DEBUG
0209 LogTrace("MuonSegFitMatrixDetails") << "[MuonSegFit::fit] p = " << p(0) << ", " << p(1) << ", " << p(2) << ", "
0210 << p(3);
0211 #endif
0212
0213
0214
0215 intercept_ = LocalPoint(p(0), p(1), 0.);
0216
0217
0218 uslope_ = p(2);
0219 vslope_ = p(3);
0220 setOutFromIP();
0221
0222
0223 setChi2();
0224
0225
0226 fitdone_ = true;
0227 }
0228
0229 void MuonSegFit::setChi2(void) {
0230 double chsq = 0.;
0231
0232 MuonRecHitContainer::const_iterator ih;
0233
0234 for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0235 LocalPoint lp = (*ih)->localPosition();
0236
0237 double u = lp.x();
0238 double v = lp.y();
0239 double z = lp.z();
0240
0241 double du = intercept_.x() + uslope_ * z - u;
0242 double dv = intercept_.y() + vslope_ * z - v;
0243
0244 SMatrixSym2 IC;
0245
0246 IC(0, 0) = (*ih)->localPositionError().xx();
0247
0248 IC(1, 0) = (*ih)->localPositionError().xy();
0249 IC(1, 1) = (*ih)->localPositionError().yy();
0250
0251
0252 #ifdef EDM_ML_DEBUG
0253 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::setChi2] IC before = \n" << IC;
0254 #endif
0255
0256
0257 bool ok = IC.Invert();
0258 if (!ok) {
0259 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::setChi2] Failed to invert covariance matrix: \n"
0260 << IC;
0261
0262 }
0263 chsq += du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
0264
0265 #ifdef EDM_ML_DEBUG
0266 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::setChi2] IC after = \n" << IC;
0267 edm::LogVerbatim("MuonSegFit")
0268 << "[GEM RecHit ] Contribution to Chi^2 of this hit :: du^2*Cov(0,0) + 2*du*dv*Cov(0,1) + dv^2*IC(1,1) = "
0269 << du * du << "*" << IC(0, 0) << " + 2.*" << du << "*" << dv << "*" << IC(0, 1) << " + " << dv * dv << "*"
0270 << IC(1, 1) << " = " << chsq;
0271 #endif
0272 }
0273
0274
0275 chi2_ = chsq;
0276 ndof_ = 2. * hits_.size() - 4;
0277
0278 #ifdef EDM_ML_DEBUG
0279 edm::LogVerbatim("MuonSegFit") << "[MuonSegFit::setChi2] chi2/ndof = " << chi2_ << "/" << ndof_;
0280 edm::LogVerbatim("MuonSegFit") << "-----------------------------------------------------";
0281
0282
0283
0284
0285 edm::LogVerbatim("MuonSegFit") << "[GEM Segment with Local Direction = " << localdir_ << " and Local Position "
0286 << intercept_ << "] can be written as:";
0287 edm::LogVerbatim("MuonSegFit") << "[ x ] = " << localdir_.x() << " * t + " << intercept_.x();
0288 edm::LogVerbatim("MuonSegFit") << "[ y ] = " << localdir_.y() << " * t + " << intercept_.y();
0289 edm::LogVerbatim("MuonSegFit") << "[ z ] = " << localdir_.z() << " * t + " << intercept_.z();
0290 edm::LogVerbatim("MuonSegFit")
0291 << "Now extrapolate to each of the GEMRecHits XY plane (so constant z = RH LP.z()) to obtain [x1,y1]";
0292 #endif
0293 }
0294
0295 MuonSegFit::SMatrixSym12 MuonSegFit::weightMatrix() {
0296 bool ok = true;
0297
0298 SMatrixSym12 matrix = ROOT::Math::SMatrixIdentity();
0299
0300 int row = 0;
0301
0302 for (MuonRecHitContainer::const_iterator it = hits_.begin(); it != hits_.end(); ++it) {
0303
0304 matrix(row, row) = scaleXError() * (*it)->localPositionError().xx();
0305 matrix(row, row + 1) = (*it)->localPositionError().xy();
0306 ++row;
0307 matrix(row, row - 1) = (*it)->localPositionError().xy();
0308 matrix(row, row) = (*it)->localPositionError().yy();
0309 ++row;
0310 }
0311
0312 ok = matrix.Invert();
0313 if (!ok) {
0314 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::weightMatrix] Failed to invert matrix: \n" << matrix;
0315
0316 }
0317 return matrix;
0318 }
0319
0320 MuonSegFit::SMatrix12by4 MuonSegFit::derivativeMatrix() {
0321 SMatrix12by4 matrix;
0322 int row = 0;
0323
0324 for (MuonRecHitContainer::const_iterator it = hits_.begin(); it != hits_.end(); ++it) {
0325 LocalPoint lp = (*it)->localPosition();
0326
0327 float z = lp.z();
0328
0329 matrix(row, 0) = 1.;
0330 matrix(row, 2) = z;
0331 ++row;
0332 matrix(row, 1) = 1.;
0333 matrix(row, 3) = z;
0334 ++row;
0335 }
0336 return matrix;
0337 }
0338
0339 void MuonSegFit::setOutFromIP() {
0340
0341
0342
0343 double dxdz = uslope_;
0344 double dydz = vslope_;
0345 double dz = 1. / sqrt(1. + dxdz * dxdz + dydz * dydz);
0346 double dx = dz * dxdz;
0347 double dy = dz * dydz;
0348 LocalVector localDir(dx, dy, dz);
0349
0350 localdir_ = (localDir).unit();
0351 #ifdef EDM_ML_DEBUG
0352 edm::LogVerbatim("MuonSegFit") << "[MuonSegFit::setOutFromIP] :: dxdz = uslope_ = " << std::setw(9) << uslope_
0353 << " dydz = vslope_ = " << std::setw(9) << vslope_ << " local dir = " << localDir;
0354 edm::LogVerbatim("MuonSegFit") << "[MuonSegFit::setOutFromIP] :: ==> local dir = " << localdir_
0355 << " localdir.phi = " << localdir_.phi();
0356 #endif
0357 }
0358
0359 AlgebraicSymMatrix MuonSegFit::covarianceMatrix() {
0360 SMatrixSym12 weights = weightMatrix();
0361 SMatrix12by4 A = derivativeMatrix();
0362 #ifdef EDM_ML_DEBUG
0363 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::covarianceMatrix] weights matrix W: \n" << weights;
0364 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::covarianceMatrix] derivatives matrix A: \n" << A;
0365 #endif
0366
0367
0368
0369
0370 bool ok;
0371 SMatrixSym4 result = ROOT::Math::SimilarityT(A, weights);
0372 #ifdef EDM_ML_DEBUG
0373 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::covarianceMatrix] (AT W A): \n" << result;
0374 #endif
0375
0376 ok = result.Invert();
0377 if (!ok) {
0378 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::calculateError] Failed to invert matrix: \n" << result;
0379
0380 }
0381 #ifdef EDM_ML_DEBUG
0382 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::covarianceMatrix] (AT W A)^-1: \n" << result;
0383 #endif
0384
0385
0386
0387 AlgebraicSymMatrix flipped = flipErrors(result);
0388
0389 return flipped;
0390 }
0391
0392 AlgebraicSymMatrix MuonSegFit::flipErrors(const SMatrixSym4& a) {
0393
0394
0395
0396
0397 #ifdef EDM_ML_DEBUG
0398 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::flipErrors] input: \n" << a;
0399 #endif
0400
0401 AlgebraicSymMatrix hold(4, 0.);
0402
0403 for (short j = 0; j != 4; ++j) {
0404 for (short i = 0; i != 4; ++i) {
0405 hold(i + 1, j + 1) = a(i, j);
0406 }
0407 }
0408
0409 #ifdef EDM_ML_DEBUG
0410 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::flipErrors] after copy:";
0411 edm::LogVerbatim("MuonSegFitMatrixDetails")
0412 << "(" << hold(1, 1) << " " << hold(1, 2) << " " << hold(1, 3) << " " << hold(1, 4);
0413 edm::LogVerbatim("MuonSegFitMatrixDetails")
0414 << " " << hold(2, 1) << " " << hold(2, 2) << " " << hold(2, 3) << " " << hold(2, 4);
0415 edm::LogVerbatim("MuonSegFitMatrixDetails")
0416 << " " << hold(3, 1) << " " << hold(3, 2) << " " << hold(3, 3) << " " << hold(3, 4);
0417 edm::LogVerbatim("MuonSegFitMatrixDetails")
0418 << " " << hold(4, 1) << " " << hold(4, 2) << " " << hold(4, 3) << " " << hold(4, 4) << ")";
0419 #endif
0420
0421
0422 hold(1, 1) = a(2, 2);
0423 hold(1, 2) = a(2, 3);
0424 hold(2, 1) = a(3, 2);
0425 hold(2, 2) = a(3, 3);
0426
0427
0428 hold(3, 3) = a(0, 0);
0429 hold(3, 4) = a(0, 1);
0430 hold(4, 3) = a(1, 0);
0431 hold(4, 4) = a(1, 1);
0432
0433
0434 hold(4, 1) = a(1, 2);
0435 hold(3, 2) = a(0, 3);
0436 hold(2, 3) = a(3, 0);
0437 hold(1, 4) = a(2, 1);
0438
0439 #ifdef EDM_ML_DEBUG
0440 edm::LogVerbatim("MuonSegFitMatrixDetails") << "[MuonSegFit::flipErrors] after flip:";
0441 edm::LogVerbatim("MuonSegFitMatrixDetails")
0442 << "(" << hold(1, 1) << " " << hold(1, 2) << " " << hold(1, 3) << " " << hold(1, 4);
0443 edm::LogVerbatim("MuonSegFitMatrixDetails")
0444 << " " << hold(2, 1) << " " << hold(2, 2) << " " << hold(2, 3) << " " << hold(2, 4);
0445 edm::LogVerbatim("MuonSegFitMatrixDetails")
0446 << " " << hold(3, 1) << " " << hold(3, 2) << " " << hold(3, 3) << " " << hold(3, 4);
0447 edm::LogVerbatim("MuonSegFitMatrixDetails")
0448 << " " << hold(4, 1) << " " << hold(4, 2) << " " << hold(4, 3) << " " << hold(4, 4) << ")";
0449 #endif
0450
0451 return hold;
0452 }
0453
0454 float MuonSegFit::xfit(float z) const {
0455
0456
0457 return intercept_.x() + uslope_ * z;
0458 }
0459
0460 float MuonSegFit::yfit(float z) const { return intercept_.y() + vslope_ * z; }
0461
0462 float MuonSegFit::xdev(float x, float z) const { return intercept_.x() + uslope_ * z - x; }
0463
0464 float MuonSegFit::ydev(float y, float z) const { return intercept_.y() + vslope_ * z - y; }
0465
0466 float MuonSegFit::Rdev(float x, float y, float z) const {
0467 return sqrt(xdev(x, z) * xdev(x, z) + ydev(y, z) * ydev(y, z));
0468 }