PixelLowPtUtilities/src/TrackFitter.cc

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 }  // namespace
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   // pt
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   // tip
0076   float valTip = charge * (center.mag() - 1 / curvature);
0077   // zip
0078   float valZip = getZip(valTip, curvature, points[0], points[1]);
0079   // phi
0080   float valPhi = getPhi(center.x(), center.y(), charge);
0081   // cot(theta), update zip
0082   float valCotTheta = getCotThetaAndUpdateZip(points[0], points[1], 1 / curvature, valPhi, valTip, valZip);
0083
0084   // errors
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   // build pixel track
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   // Recalculate ZIP
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   // phi = asin(r*rho/2) with asin(x) ~= x+x**3/(2*3)
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   {                       // transverse
0161     float c_ms = 0.0115;  //0.0115;
0162     float s_le = 0.0095;  //0.0123;
0163     float s_ms2 = c_ms * c_ms / (pt * pt) * coshEta;
0164
0165     errTip = sqrt(s_le * s_le + s_ms2);
0166   }
0167
0168   {  // z
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 }