File indexing completed on 2023-03-31 23:02:12
0001 #include "CommonTools/Statistics/interface/LinearFit.h"
0002 #include "CommonTools/Utils/interface/DynArray.h"
0003 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0004 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0005 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0006 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0007 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0011 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0012 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0013 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0014 #include "MagneticField/Engine/interface/MagneticField.h"
0015 #include "RecoTracker/PixelLowPtUtilities/interface/TrackFitter.h"
0016 #include "RecoTracker/PixelTrackFitting/interface/CircleFromThreePoints.h"
0017 #include "RecoTracker/PixelTrackFitting/interface/PixelTrackBuilder.h"
0018 #include "RecoTracker/PixelTrackFitting/interface/RZLine.h"
0019 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0020 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0021 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0022
0023 using namespace std;
0024
0025 namespace {
0026
0027 int getCharge(const DynArray<GlobalPoint>& points) {
0028 GlobalVector v21 = points[1] - points[0];
0029 GlobalVector v32 = points[2] - points[1];
0030 float dphi = v32.phi() - v21.phi();
0031 while (dphi > M_PI)
0032 dphi -= 2 * M_PI;
0033 while (dphi < -M_PI)
0034 dphi += 2 * M_PI;
0035 return (dphi > 0) ? -1 : 1;
0036 }
0037
0038 }
0039
0040
0041 std::unique_ptr<reco::Track> TrackFitter::run(const std::vector<const TrackingRecHit*>& hits,
0042 const TrackingRegion& region) const {
0043 std::unique_ptr<reco::Track> ret;
0044
0045 int nhits = hits.size();
0046 if (nhits < 2)
0047 return ret;
0048
0049 declareDynArray(GlobalPoint, nhits, points);
0050 declareDynArray(GlobalError, nhits, errors);
0051 declareDynArray(bool, nhits, isBarrel);
0052
0053 unsigned int i = 0;
0054 for (auto const& ih : hits) {
0055 auto recHit = theTTRecHitBuilder->build(ih);
0056
0057 points[i] = recHit->globalPosition();
0058 errors[i] = recHit->globalPositionError();
0059 isBarrel[++i] = recHit->detUnit()->type().isBarrel();
0060 }
0061
0062 CircleFromThreePoints circle = (nhits == 2) ? CircleFromThreePoints(GlobalPoint(0., 0., 0.), points[0], points[1])
0063 : CircleFromThreePoints(points[0], points[1], points[2]);
0064
0065 int charge = getCharge(points);
0066 float curvature = circle.curvature();
0067
0068
0069 float invPt = PixelRecoUtilities::inversePt(curvature, *theField);
0070 float valPt = (invPt > 1.e-4) ? 1. / invPt : 1.e4;
0071 float errPt = 0.055 * valPt + 0.017 * valPt * valPt;
0072
0073 CircleFromThreePoints::Vector2D center = circle.center();
0074
0075
0076 float valTip = charge * (center.mag() - 1 / curvature);
0077
0078 float valZip = getZip(valTip, curvature, points[0], points[1]);
0079
0080 float valPhi = getPhi(center.x(), center.y(), charge);
0081
0082 float valCotTheta = getCotThetaAndUpdateZip(points[0], points[1], 1 / curvature, valPhi, valTip, valZip);
0083
0084
0085 float errTip, errZip;
0086 getErrTipAndErrZip(valPt, points.back().eta(), errTip, errZip);
0087 float errPhi = 0.002;
0088 float errCotTheta = 0.002;
0089
0090 float chi2 = 0;
0091 if (nhits > 2) {
0092 RZLine rzLine(points, errors, isBarrel);
0093 chi2 = rzLine.chi2();
0094 }
0095
0096
0097 PixelTrackBuilder builder;
0098
0099 Measurement1D pt(valPt, errPt);
0100 Measurement1D phi(valPhi, errPhi);
0101 Measurement1D cotTheta(valCotTheta, errCotTheta);
0102 Measurement1D tip(valTip, errTip);
0103 Measurement1D zip(valZip, errZip);
0104
0105 ret.reset(builder.build(pt, phi, cotTheta, tip, zip, chi2, charge, hits, theField));
0106 return ret;
0107 }
0108
0109
0110 float TrackFitter::getCotThetaAndUpdateZip(
0111 const GlobalPoint& inner, const GlobalPoint& outer, float radius, float phi, float d0, float& zip) const {
0112 float chi = phi - M_PI_2;
0113 GlobalPoint IP(d0 * cos(chi), d0 * sin(chi), zip);
0114
0115 float phi1 = 2 * asin(0.5 * (inner - IP).perp() / radius);
0116 float phi2 = 2 * asin(0.5 * (outer - IP).perp() / radius);
0117
0118 float dr = radius * (phi2 - phi1);
0119 float dz = outer.z() - inner.z();
0120
0121
0122 zip = (inner.z() * phi2 - outer.z() * phi1) / (phi2 - phi1);
0123
0124 return (fabs(dr) > 1.e-3) ? dz / dr : 0;
0125 }
0126
0127
0128 float TrackFitter::getPhi(float xC, float yC, int charge) const {
0129 float phiC;
0130
0131 if (charge > 0)
0132 phiC = atan2(xC, -yC);
0133 else
0134 phiC = atan2(-xC, yC);
0135
0136 return phiC;
0137 }
0138
0139
0140 float TrackFitter::getZip(float d0, float curv, const GlobalPoint& inner, const GlobalPoint& outer) const {
0141
0142 float rho3 = curv * curv * curv;
0143
0144 float r1 = inner.perp();
0145 double phi1 = r1 * curv / 2 + inner.perp2() * r1 * rho3 / 48.;
0146
0147 float r2 = outer.perp();
0148 double phi2 = r2 * curv / 2 + outer.perp2() * r2 * rho3 / 48.;
0149
0150 double z1 = inner.z();
0151 double z2 = outer.z();
0152
0153 return z1 - phi1 / (phi1 - phi2) * (z1 - z2);
0154 }
0155
0156
0157 void TrackFitter::getErrTipAndErrZip(float pt, float eta, float& errTip, float& errZip) const {
0158 float coshEta = cosh(eta);
0159
0160 {
0161 float c_ms = 0.0115;
0162 float s_le = 0.0095;
0163 float s_ms2 = c_ms * c_ms / (pt * pt) * coshEta;
0164
0165 errTip = sqrt(s_le * s_le + s_ms2);
0166 }
0167
0168 {
0169 float c_ms = 0.0070;
0170 float s_le = 0.0135;
0171
0172 errZip = sqrt((s_le * s_le + c_ms * c_ms / (pt * pt)) * coshEta * coshEta * coshEta);
0173 }
0174 }