File indexing completed on 2024-04-06 12:25:59
0001
0002
0003
0004 #include "RecoLocalMuon/CSCSegment/src/CSCCondSegFit.h"
0005
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 #include <algorithm>
0011
0012 void CSCCondSegFit::fit(bool condpass1, bool condpass2) {
0013
0014
0015
0016
0017 lex_.clear();
0018
0019 if (condpass1 && !condpass2) {
0020 correctTheCovX();
0021 if (lex_.size() != hits_.size()) {
0022 LogTrace("CSCSegment|CSC") << "[CSCSegFit::fit] lex_.size()!=hits_.size() ALARM! Please inform CSC DPG "
0023 << std::endl;
0024 }
0025 }
0026
0027 SMatrix4 M;
0028 SVector4 B;
0029
0030 CSCSetOfHits::const_iterator ih = hits_.begin();
0031
0032 for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0033 const CSCRecHit2D& hit = (**ih);
0034 const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer());
0035 GlobalPoint gp = layer->toGlobal(hit.localPosition());
0036 LocalPoint lp = chamber()->toLocal(gp);
0037
0038
0039 double u = lp.x();
0040 double v = lp.y();
0041 double z = lp.z();
0042
0043
0044 SMatrixSym2 IC;
0045
0046
0047 if (condpass1 && !condpass2) {
0048 IC(0, 0) = lex_.at(ih - hits_.begin());
0049 } else {
0050 IC(0, 0) = hit.localPositionError().xx();
0051 }
0052 IC(1, 1) = hit.localPositionError().yy();
0053
0054
0055 IC(1, 0) = hit.localPositionError().xy();
0056
0057
0058
0059 if (condpass2) {
0060 correctTheCovMatrix(IC);
0061 }
0062
0063
0064 bool ok = IC.Invert();
0065 if (!ok) {
0066 LogTrace("CSCSegment|CSC") << "[CSCSegFit::fit] Failed to invert covariance matrix: \n" << IC;
0067
0068 }
0069
0070 M(0, 0) += IC(0, 0);
0071 M(0, 1) += IC(0, 1);
0072 M(0, 2) += IC(0, 0) * z;
0073 M(0, 3) += IC(0, 1) * z;
0074 B(0) += u * IC(0, 0) + v * IC(0, 1);
0075
0076 M(1, 0) += IC(1, 0);
0077 M(1, 1) += IC(1, 1);
0078 M(1, 2) += IC(1, 0) * z;
0079 M(1, 3) += IC(1, 1) * z;
0080 B(1) += u * IC(1, 0) + v * IC(1, 1);
0081
0082 M(2, 0) += IC(0, 0) * z;
0083 M(2, 1) += IC(0, 1) * z;
0084 M(2, 2) += IC(0, 0) * z * z;
0085 M(2, 3) += IC(0, 1) * z * z;
0086 B(2) += (u * IC(0, 0) + v * IC(0, 1)) * z;
0087
0088 M(3, 0) += IC(1, 0) * z;
0089 M(3, 1) += IC(1, 1) * z;
0090 M(3, 2) += IC(1, 0) * z * z;
0091 M(3, 3) += IC(1, 1) * z * z;
0092 B(3) += (u * IC(1, 0) + v * IC(1, 1)) * z;
0093 }
0094
0095 SVector4 p;
0096 bool ok = M.Invert();
0097 if (!ok) {
0098 LogTrace("CSCSegment|CSC") << "[CSCSegFit::fit] Failed to invert matrix: \n" << M;
0099
0100 } else {
0101 p = M * B;
0102 }
0103
0104
0105
0106
0107
0108 intercept_ = LocalPoint(p(0), p(1), 0.);
0109
0110
0111 uslope_ = p(2);
0112 vslope_ = p(3);
0113 setOutFromIP();
0114
0115
0116 setChi2(condpass1, condpass2);
0117 }
0118
0119 void CSCCondSegFit::setChi2(bool condpass1, bool condpass2) {
0120 double chsq = 0.;
0121
0122 CSCSetOfHits::const_iterator ih;
0123 for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0124 const CSCRecHit2D& hit = (**ih);
0125 const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer());
0126 GlobalPoint gp = layer->toGlobal(hit.localPosition());
0127 LocalPoint lp = chamber()->toLocal(gp);
0128
0129 double u = lp.x();
0130 double v = lp.y();
0131 double z = lp.z();
0132
0133 double du = intercept_.x() + uslope_ * z - u;
0134 double dv = intercept_.y() + vslope_ * z - v;
0135
0136
0137
0138 SMatrixSym2 IC;
0139
0140 if (condpass1 && !condpass2) {
0141 IC(0, 0) = lex_.at(ih - hits_.begin());
0142 } else {
0143 IC(0, 0) = hit.localPositionError().xx();
0144 }
0145
0146 IC(1, 0) = hit.localPositionError().xy();
0147 IC(1, 1) = hit.localPositionError().yy();
0148
0149
0150
0151
0152
0153 if (condpass2) {
0154 correctTheCovMatrix(IC);
0155 }
0156
0157
0158 bool ok = IC.Invert();
0159 if (!ok) {
0160 LogTrace("CSCSegment|CSC") << "[CSCSegFit::setChi2] Failed to invert covariance matrix: \n" << IC;
0161
0162 }
0163
0164 chsq += du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
0165 }
0166
0167
0168 chi2_ = chsq;
0169 ndof_ = 2. * hits_.size() - 4;
0170
0171
0172 }
0173
0174 void CSCCondSegFit::correctTheCovX(void) {
0175 std::vector<double> uu, vv, zz;
0176
0177 lex_.clear();
0178
0179 double sum_U_err = 0.0;
0180 double sum_Z_U_err = 0.0;
0181 double sum_Z2_U_err = 0.0;
0182 double sum_U_U_err = 0.0;
0183 double sum_UZ_U_err = 0.0;
0184 std::vector<double> chiUZind;
0185 std::vector<double>::iterator chiContribution;
0186 double chiUZ = 0.0;
0187 CSCSetOfHits::const_iterator ih = hits_.begin();
0188 for (ih = hits_.begin(); ih != hits_.end(); ++ih) {
0189 const CSCRecHit2D& hit = (**ih);
0190 lex_.push_back(hit.localPositionError().xx());
0191
0192 const CSCLayer* layer = chamber()->layer(hit.cscDetId().layer());
0193 GlobalPoint gp = layer->toGlobal(hit.localPosition());
0194 LocalPoint lp = chamber()->toLocal(gp);
0195
0196 double u = lp.x();
0197 double v = lp.y();
0198 double z = lp.z();
0199 uu.push_back(u);
0200 vv.push_back(v);
0201 zz.push_back(z);
0202
0203 sum_U_err += 1. / lex_.back();
0204 sum_Z_U_err += z / lex_.back();
0205 sum_Z2_U_err += (z * z) / lex_.back();
0206 sum_U_U_err += u / lex_.back();
0207 sum_UZ_U_err += (u * z) / lex_.back();
0208 }
0209
0210
0211
0212
0213 double denom = sum_U_err * sum_Z2_U_err - pow(sum_Z_U_err, 2);
0214 double U0 = (sum_Z2_U_err * sum_U_U_err - sum_Z_U_err * sum_UZ_U_err) / denom;
0215 double UZ = (sum_U_err * sum_UZ_U_err - sum_Z_U_err * sum_U_U_err) / denom;
0216
0217
0218
0219
0220 for (unsigned i = 0; i < uu.size(); ++i) {
0221 double uMean = U0 + UZ * zz[i];
0222 chiUZind.push_back((pow((uMean - uu[i]), 2)) / lex_[i]);
0223 chiUZ += (pow((uMean - uu[i]), 2)) / lex_[i];
0224 }
0225 chiUZ = chiUZ / (uu.size() - 2);
0226
0227 if (chiUZ >= chi2Norm_) {
0228 double chi2uCorrection = chiUZ / chi2Norm_;
0229 for (unsigned i = 0; i < uu.size(); ++i)
0230 lex_[i] = lex_[i] * chi2uCorrection;
0231 setScaleXError(chi2uCorrection);
0232 }
0233
0234
0235 chiContribution = max_element(chiUZind.begin(), chiUZind.end());
0236 worstHit_ = chiContribution - chiUZind.begin();
0237 }
0238
0239 void CSCCondSegFit::correctTheCovMatrix(SMatrixSym2& IC) {
0240
0241 double condNumberCorr2 = 0.0;
0242 double detCov = 0.0;
0243 double diag1 = 0.0;
0244 double diag2 = 0.0;
0245 double IC_01_corr = 0.0;
0246 double IC_00_corr = 0.0;
0247 if (!covToAnyNumberAll_) {
0248
0249 condNumberCorr2 = condSeed2_ * IC(1, 1);
0250 diag1 = IC(0, 0) * IC(1, 1);
0251 diag2 = IC(0, 1) * IC(0, 1);
0252 detCov = fabs(diag1 - diag2);
0253 if ((diag1 < condNumberCorr2) && (diag2 < condNumberCorr2)) {
0254 if (covToAnyNumber_)
0255 IC(0, 1) = covAnyNumber_;
0256 else {
0257 IC_00_corr = condSeed1_ + fabs(IC(0, 1)) / IC(1, 1);
0258 IC(0, 0) = IC_00_corr;
0259 }
0260 }
0261
0262 if (((detCov < condNumberCorr2) && (diag1 > condNumberCorr2)) ||
0263 ((diag2 > condNumberCorr2) && (detCov < condNumberCorr2))) {
0264 if (covToAnyNumber_)
0265 IC(0, 1) = covAnyNumber_;
0266 else {
0267 IC_01_corr = sqrt(fabs(diag1 - condNumberCorr2));
0268 if (IC(0, 1) < 0)
0269 IC(0, 1) = (-1) * IC_01_corr;
0270 else
0271 IC(0, 1) = IC_01_corr;
0272 }
0273 }
0274 } else {
0275 IC(0, 1) = covAnyNumber_;
0276 }
0277
0278
0279
0280 IC(1, 0) = IC(0, 1);
0281 }