Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:30:42

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     double chiSq = 0.0;
0112 
0113     unsigned int j = 0;
0114 
0115     largestresid_ = -1.0;
0116     ilargestresid_ = -1;
0117 
0118     // Scattering term scaling as 1/Pt.
0119     double sigmaScat = settings_->kalmanMultiScattTerm() * std::abs(qOverPt_seed_);
0120 
0121     for (unsigned int i = 0; i < n; i++) {
0122       double ri = stubs_[i]->r();
0123       double zi = stubs_[i]->z();
0124       double phii = stubs_[i]->phi();
0125       double sigmaPerp = stubs_[i]->sigmaPerp();
0126       double sigmaPar = stubs_[i]->sigmaPar();
0127       sigmaPerp = sqrt(sigmaPerp * sigmaPerp + sigmaScat * sigmaScat * ri * ri);
0128 
0129       if (stubs_[i]->barrel()) {
0130         double halfRinvRi = 0.5 * ri * rInv;
0131         double aSinHalfRinvRi = halfRinvRi + (2. / 3.) * halfRinvRi * halfRinvRi * halfRinvRi;
0132         double deltaphi = reco::deltaPhi(phi0 - aSinHalfRinvRi - phii, 0.);
0133         delta[j++] = (ri * deltaphi) / sigmaPerp;
0134         delta[j++] = (z0 + (2.0 / rInv) * t * aSinHalfRinvRi - zi) / sigmaPar;
0135       } else {
0136         double tInv = 1 / t;
0137         double r_track = (zi - z0) * tInv;
0138         double phi_track = phi0 - 0.5 * rInv * (zi - z0) * tInv;
0139         int iphi = stubs_[i]->iphi();
0140 
0141         // N.B. These represent HALF the width and number of strips of sensor.
0142         double width = 0.5 * stubs_[i]->trackerModule()->sensorWidth();
0143         double nstrip = 0.5 * stubs_[i]->nStrips();
0144 
0145         double Deltai = width * (iphi - nstrip) / nstrip;  // Non-radial endcap 2S strip correction
0146 
0147         if (stubs_[i]->z() > 0.0)
0148           Deltai = -Deltai;
0149 
0150         double DeltaiOverRi = Deltai / ri;
0151         double theta0 = DeltaiOverRi + (2. / 3.) * DeltaiOverRi * DeltaiOverRi * DeltaiOverRi;
0152         double Delta = Deltai - r_track * (theta0 - (phi_track - phii));
0153 
0154         delta[j++] = (r_track - ri) / sigmaPar;
0155         delta[j++] = Delta / sigmaPerp;
0156       }
0157 
0158       chiSq += delta[j - 2] * delta[j - 2] + delta[j - 1] * delta[j - 1];
0159 
0160       if (std::abs(delta[j - 2]) > largestresid_) {
0161         largestresid_ = std::abs(delta[j - 2]);
0162         ilargestresid_ = i;
0163       }
0164 
0165       if (std::abs(delta[j - 1]) > largestresid_) {
0166         largestresid_ = std::abs(delta[j - 1]);
0167         ilargestresid_ = i;
0168       }
0169     }
0170 
0171     return delta;
0172   }
0173 
0174 }  // namespace tmtt