PixelTrackFitting/src/PixelFitterByHelixProjections.cc

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 //#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
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     // the cross product will tell me...
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       GlobalVector v21 = points[1]-points[0];
0047       GlobalVector v32 = points[2]-points[1];
0048       float dphi = v32.phi() - v21.phi();
0049       while (dphi >  Geom::fpi()) dphi -=  Geom::ftwoPi();
0050       while (dphi < -Geom::fpi()) dphi +=  Geom::ftwoPi();
0051       return (dphi > 0) ? -1 : 1;
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     //phi = asin(r*rho/2) with asin(x) ~= x+x**3/(2*3) = x(1+x*x/6);
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 }  // namespace
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   // Rescale down the error to take into accont the fact that the
0141   // inner pixel barrel layer for PhaseI is closer to the interaction
0142   // point. The effective scale factor has been derived by checking
0143   // that the pulls of the pixelVertices derived from the pixelTracks
0144   // have the correct mean and sigma.
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 }