File indexing completed on 2023-03-31 23:02:23
0001 #include "RecoTracker/PixelTrackFitting/interface/PixelFitterByHelixProjections.h"
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005
0006 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0007 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0008 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0009
0010 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0011 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0012 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0013
0014 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0015
0016 #include "CommonTools/Statistics/interface/LinearFit.h"
0017 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0018
0019 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0020 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0021
0022 #include "MagneticField/Engine/interface/MagneticField.h"
0023
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025
0026 #include "RecoTracker/PixelTrackFitting/interface/RZLine.h"
0027 #include "RecoTracker/PixelTrackFitting/interface/CircleFromThreePoints.h"
0028 #include "RecoTracker/PixelTrackFitting/interface/PixelTrackBuilder.h"
0029 #include "RecoTracker/PixelTrackFitting/interface/PixelTrackErrorParam.h"
0030 #include "DataFormats/GeometryVector/interface/Pi.h"
0031
0032 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0033
0034 #include "CommonTools/Utils/interface/DynArray.h"
0035
0036 using namespace std;
0037
0038 namespace {
0039
0040 int charge(DynArray<GlobalPoint> const& points) {
0041
0042 float dir = (points[1].x() - points[0].x()) * (points[2].y() - points[1].y()) -
0043 (points[1].y() - points[0].y()) * (points[2].x() - points[1].x());
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053 return (dir > 0) ? -1 : 1;
0054 }
0055
0056 float cotTheta(const GlobalPoint& inner, const GlobalPoint& outer) {
0057 float dr = outer.perp() - inner.perp();
0058 float dz = outer.z() - inner.z();
0059 return (std::abs(dr) > 1.e-3f) ? dz / dr : 0;
0060 }
0061
0062 inline float phi(float xC, float yC, int charge) { return (charge > 0) ? std::atan2(xC, -yC) : std::atan2(-xC, yC); }
0063
0064 float zip(float d0, float phi_p, float curv, const GlobalPoint& pinner, const GlobalPoint& pouter) {
0065
0066
0067
0068
0069 float phi0 = phi_p - Geom::fhalfPi();
0070 GlobalPoint pca(d0 * std::cos(phi0), d0 * std::sin(phi0), 0.);
0071
0072 constexpr float o24 = 1.f / 24.f;
0073 float rho2 = curv * curv;
0074 float r1s = (pinner - pca).perp2();
0075 double phi1 = std::sqrt(r1s) * (curv * 0.5f) * (1.f + r1s * (rho2 * o24));
0076 float r2s = (pouter - pca).perp2();
0077 double phi2 = std::sqrt(r2s) * (curv * 0.5f) * (1.f + r2s * (rho2 * o24));
0078 double z1 = pinner.z();
0079 double z2 = pouter.z();
0080
0081 if (fabs(curv) > 1.e-5)
0082 return z1 - phi1 / (phi1 - phi2) * (z1 - z2);
0083 else {
0084 double dr = std::max(std::sqrt(r2s) - std::sqrt(r1s), 1.e-5f);
0085 return z1 - std::sqrt(r1s) * (z2 - z1) / dr;
0086 }
0087 }
0088 }
0089
0090 PixelFitterByHelixProjections::PixelFitterByHelixProjections(const TrackerTopology* ttopo,
0091 const MagneticField* field,
0092 bool scaleErrorsForBPix1,
0093 float scaleFactor)
0094 : theTopo(ttopo), theField(field), thescaleErrorsForBPix1(scaleErrorsForBPix1), thescaleFactor(scaleFactor) {}
0095
0096 std::unique_ptr<reco::Track> PixelFitterByHelixProjections::run(const std::vector<const TrackingRecHit*>& hits,
0097 const TrackingRegion& region) const {
0098 std::unique_ptr<reco::Track> ret;
0099
0100 int nhits = hits.size();
0101 if (nhits < 2)
0102 return ret;
0103
0104 declareDynArray(GlobalPoint, nhits, points);
0105 declareDynArray(GlobalError, nhits, errors);
0106 declareDynArray(bool, nhits, isBarrel);
0107
0108 for (int i = 0; i != nhits; ++i) {
0109 auto const& recHit = hits[i];
0110 points[i] = GlobalPoint(recHit->globalPosition().basicVector() - region.origin().basicVector());
0111 errors[i] = recHit->globalPositionError();
0112 isBarrel[i] = recHit->detUnit()->type().isBarrel();
0113 }
0114
0115 CircleFromThreePoints circle = (nhits == 2) ? CircleFromThreePoints(GlobalPoint(0., 0., 0.), points[0], points[1])
0116 : CircleFromThreePoints(points[0], points[1], points[2]);
0117
0118 float valPhi, valTip, valPt;
0119
0120 int iCharge = charge(points);
0121 float curvature = circle.curvature();
0122
0123 if ((curvature > 1.e-4) && (LIKELY(theField->inverseBzAtOriginInGeV()) > 0.01)) {
0124 float invPt = PixelRecoUtilities::inversePt(circle.curvature(), *theField);
0125 valPt = (invPt > 1.e-4f) ? 1.f / invPt : 1.e4f;
0126 CircleFromThreePoints::Vector2D center = circle.center();
0127 valTip = iCharge * (center.mag() - 1.f / curvature);
0128 valPhi = phi(center.x(), center.y(), iCharge);
0129 } else {
0130 valPt = 1.e4f;
0131 GlobalVector direction(points[1] - points[0]);
0132 valPhi = direction.barePhi();
0133 valTip = -points[0].x() * sin(valPhi) + points[0].y() * cos(valPhi);
0134 }
0135
0136 float valCotTheta = cotTheta(points[0], points[1]);
0137 float valEta = std::asinh(valCotTheta);
0138 float valZip = zip(valTip, valPhi, curvature, points[0], points[1]);
0139
0140
0141
0142
0143
0144
0145 float errFactor = 1.;
0146 if (thescaleErrorsForBPix1 && (hits[0]->geographicalId().subdetId() == PixelSubdetector::PixelBarrel) &&
0147 (theTopo->pxbLayer(hits[0]->geographicalId()) == 1))
0148 errFactor = thescaleFactor;
0149
0150 PixelTrackErrorParam param(valEta, valPt);
0151 float errValPt = errFactor * param.errPt();
0152 float errValCot = errFactor * param.errCot();
0153 float errValTip = errFactor * param.errTip();
0154 float errValPhi = errFactor * param.errPhi();
0155 float errValZip = errFactor * param.errZip();
0156
0157 float chi2 = 0;
0158 if (nhits > 2) {
0159 RZLine rzLine(points, errors, isBarrel);
0160 chi2 = rzLine.chi2();
0161 }
0162
0163 PixelTrackBuilder builder;
0164 Measurement1D pt(valPt, errValPt);
0165 Measurement1D phi(valPhi, errValPhi);
0166 Measurement1D cotTheta(valCotTheta, errValCot);
0167 Measurement1D tip(valTip, errValTip);
0168 Measurement1D zip(valZip, errValZip);
0169
0170 ret.reset(builder.build(pt, phi, cotTheta, tip, zip, chi2, iCharge, hits, theField, region.origin()));
0171 return ret;
0172 }