TrackFindingTMTT/src/ChiSquaredFit4.cc

0001 #include "L1Trigger/TrackFindingTMTT/interface/ChiSquaredFit4.h"
0002 #include "DataFormats/Math/interface/deltaPhi.h"
0003
0004 using namespace std;
0005
0006 namespace tmtt {
0007
0008   ChiSquaredFit4::ChiSquaredFit4(const Settings* settings, const uint nPar) : ChiSquaredFitBase(settings, nPar) {
0009     largestresid_ = -1.0;
0010     ilargestresid_ = -1;
0011   }
0012
0013   TVectorD ChiSquaredFit4::seed(const L1track3D& l1track3D) {
0014     TVectorD x(4);
0015     x[INVR] = settings_->invPtToInvR() * l1track3D.qOverPt();
0016     x[PHI0] = l1track3D.phi0();
0017     x[T] = l1track3D.tanLambda();
0018     x[Z0] = l1track3D.z0();
0019     return x;
0020   }
0021
0022   //=== Calculate derivatives of track intercept with respect to track parameters
0023
0024   TMatrixD ChiSquaredFit4::D(const TVectorD& x) {
0025     TMatrixD D(2 * stubs_.size(), nPar_);  // Empty matrix
0026     D.Zero();
0027     int j = 0;
0028     double rInv = x[INVR];
0029     double phi0 = x[PHI0];
0030     double t = x[T];
0031     for (unsigned i = 0; i < stubs_.size(); i++) {
0032       double ri = stubs_[i]->r();
0033       if (stubs_[i]->barrel()) {
0034         // Derivatives of r*phi
0035         D(j, INVR) = -0.5 * ri * ri;
0036         D(j, PHI0) = ri;
0037         j++;
0038         // Derivatives of z
0039         D(j, T) = ri;
0040         D(j, Z0) = 1;
0041         j++;
0042       } else {
0043         double phii = stubs_[i]->phi();
0044         int iphi = stubs_[i]->iphi();
0045
0046         // N.B. These represent HALF the width and number of strips of sensor.
0047         double width = 0.5 * stubs_[i]->trackerModule()->sensorWidth();
0048         double nstrip = 0.5 * stubs_[i]->nStrips();
0049
0050         double Deltai = width * (iphi - nstrip) / nstrip;  // Non-radial endcap 2S strip correction
0051         if (stubs_[i]->z() > 0.0)
0052           Deltai = -Deltai;
0053         double DeltaiOverRi = Deltai / ri;
0054         double theta0 = DeltaiOverRi + (2. / 3.) * DeltaiOverRi * DeltaiOverRi * DeltaiOverRi;
0055
0056         double phi_track = phi0 - 0.5 * rInv * ri;  //Expected phi hit given the track
0057
0058         double tInv = 1 / t;
0059         // Derivatives of r
0060         D(j, INVR) = -1 * ri * ri * ri * rInv;
0061         D(j, PHI0) = 0;
0062         D(j, T) = -ri * tInv;
0063         D(j, Z0) = -1 * tInv;
0064         j++;
0065         // Derivatives of r*phi
0066         D(j, INVR) = -0.5 * ri * ri;
0067         D(j, PHI0) = ri;
0068         D(j, T) = ri * 0.5 * rInv * ri * tInv - ((phi_track - phii) - theta0) * ri * tInv;
0069         D(j, Z0) = ri * 0.5 * rInv * tInv - ((phi_track - phii) - theta0) * tInv;
0070         j++;
0071       }
0072     }
0073     return D;
0074   }
0075
0076   //=== In principle, this is the stub position covariance matrix.
0077   //=== In practice, it misses a factor "sigma", because unconventionally, this is absorbed into residuals() function.
0078
0079   TMatrixD ChiSquaredFit4::Vinv() {
0080     TMatrixD Vinv(2 * stubs_.size(), 2 * stubs_.size());
0081     // Scattering term scaling as 1/Pt.
0082     double sigmaScat = settings_->kalmanMultiScattTerm() * std::abs(qOverPt_seed_);
0083     for (unsigned i = 0; i < stubs_.size(); i++) {
0084       double sigmaPerp = stubs_[i]->sigmaPerp();
0085       double sigmaPar = stubs_[i]->sigmaPar();
0086       double ri = stubs_[i]->r();
0087       sigmaPerp = sqrt(sigmaPerp * sigmaPerp + sigmaScat * sigmaScat * ri * ri);
0088       if (stubs_[i]->barrel()) {
0089         Vinv(2 * i, 2 * i) = 1 / sigmaPerp;
0090         Vinv(2 * i + 1, 2 * i + 1) = 1 / sigmaPar;
0091       } else {
0092         Vinv(2 * i, 2 * i) = 1 / sigmaPar;
0093         Vinv(2 * i + 1, 2 * i + 1) = 1 / sigmaPerp;
0094       }
0095     }
0096     return Vinv;
0097   }
0098
0099   //=== Calculate residuals w.r.t. track divided by uncertainty.
0100
0101   TVectorD ChiSquaredFit4::residuals(const TVectorD& x) {
0102     unsigned int n = stubs_.size();
0103
0104     TVectorD delta(2 * n);
0105
0106     double rInv = x[INVR];
0107     double phi0 = x[PHI0];
0108     double t = x[T];
0109     double z0 = x[Z0];
0110
0111     unsigned int j = 0;
0112
0113     largestresid_ = -1.0;
0114     ilargestresid_ = -1;
0115
0116     // Scattering term scaling as 1/Pt.
0117     double sigmaScat = settings_->kalmanMultiScattTerm() * std::abs(qOverPt_seed_);
0118
0119     for (unsigned int i = 0; i < n; i++) {
0120       double ri = stubs_[i]->r();
0121       double zi = stubs_[i]->z();
0122       double phii = stubs_[i]->phi();
0123       double sigmaPerp = stubs_[i]->sigmaPerp();
0124       double sigmaPar = stubs_[i]->sigmaPar();
0125       sigmaPerp = sqrt(sigmaPerp * sigmaPerp + sigmaScat * sigmaScat * ri * ri);
0126
0127       if (stubs_[i]->barrel()) {
0128         double halfRinvRi = 0.5 * ri * rInv;
0129         double aSinHalfRinvRi = halfRinvRi + (2. / 3.) * halfRinvRi * halfRinvRi * halfRinvRi;
0130         double deltaphi = reco::deltaPhi(phi0 - aSinHalfRinvRi - phii, 0.);
0131         delta[j++] = (ri * deltaphi) / sigmaPerp;
0132         delta[j++] = (z0 + (2.0 / rInv) * t * aSinHalfRinvRi - zi) / sigmaPar;
0133       } else {
0134         double tInv = 1 / t;
0135         double r_track = (zi - z0) * tInv;
0136         double phi_track = phi0 - 0.5 * rInv * (zi - z0) * tInv;
0137         int iphi = stubs_[i]->iphi();
0138
0139         // N.B. These represent HALF the width and number of strips of sensor.
0140         double width = 0.5 * stubs_[i]->trackerModule()->sensorWidth();
0141         double nstrip = 0.5 * stubs_[i]->nStrips();
0142
0143         double Deltai = width * (iphi - nstrip) / nstrip;  // Non-radial endcap 2S strip correction
0144
0145         if (stubs_[i]->z() > 0.0)
0146           Deltai = -Deltai;
0147
0148         double DeltaiOverRi = Deltai / ri;
0149         double theta0 = DeltaiOverRi + (2. / 3.) * DeltaiOverRi * DeltaiOverRi * DeltaiOverRi;
0150         double Delta = Deltai - r_track * (theta0 - (phi_track - phii));
0151
0152         delta[j++] = (r_track - ri) / sigmaPar;
0153         delta[j++] = Delta / sigmaPerp;
0154       }
0155
0156       if (std::abs(delta[j - 2]) > largestresid_) {
0157         largestresid_ = std::abs(delta[j - 2]);
0158         ilargestresid_ = i;
0159       }
0160
0161       if (std::abs(delta[j - 1]) > largestresid_) {
0162         largestresid_ = std::abs(delta[j - 1]);
0163         ilargestresid_ = i;
0164       }
0165     }
0166
0167     return delta;
0168   }
0169
0170 }  // namespace tmtt