PixelTrackFitting/src/KFBasedPixelFitter.cc

0001 #include "RecoTracker/PixelTrackFitting/interface/KFBasedPixelFitter.h"
0002
0003 #include "FWCore/Framework/interface/Event.h"
0004
0005 #include "MagneticField/Engine/interface/MagneticField.h"
0006
0007 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0008 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0009
0010 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0011
0012 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0013 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0015
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017
0018 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0019 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0020
0021 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0022 #include "RecoTracker/PixelTrackFitting/interface/CircleFromThreePoints.h"
0023
0024 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0025
0026 template <class T>
0027 inline T sqr(T t) {
0028   return t * t;
0029 }
0030
0031 KFBasedPixelFitter::MyBeamSpotHit::MyBeamSpotHit(const reco::BeamSpot &beamSpot, const GeomDet *geom)
0032     : TValidTrackingRecHit(*geom) {
0033   localPosition_ = LocalPoint(0., 0., 0.);
0034   //neglect XY differences and BS slope
0035   localError_ = LocalError(sqr(beamSpot.BeamWidthX()), 0.0, sqr(beamSpot.sigmaZ()));
0036 }
0037
0038 AlgebraicVector KFBasedPixelFitter::MyBeamSpotHit::parameters() const {
0039   AlgebraicVector result(1);
0040   result[0] = localPosition().x();
0041   return result;
0042 }
0043 AlgebraicSymMatrix KFBasedPixelFitter::MyBeamSpotHit::parametersError() const {
0044   LocalError le = localPositionError();
0045   AlgebraicSymMatrix m(1);
0046   m[0][0] = le.xx();
0047   return m;
0048 }
0049 AlgebraicMatrix KFBasedPixelFitter::MyBeamSpotHit::projectionMatrix() const {
0050   AlgebraicMatrix matrix(1, 5, 0);
0051   matrix[0][3] = 1;
0052   return matrix;
0053 }
0054
0055 KFBasedPixelFitter::KFBasedPixelFitter(const Propagator *propagator,
0056                                        const Propagator *opropagator,
0057                                        const TransientTrackingRecHitBuilder *ttrhBuilder,
0058                                        const TrackerGeometry *tracker,
0059                                        const MagneticField *field,
0060                                        const reco::BeamSpot *beamSpot)
0061     : thePropagator(propagator),
0062       theOPropagator(opropagator),
0063       theTTRHBuilder(ttrhBuilder),
0064       theTracker(tracker),
0065       theField(field),
0066       theBeamSpot(beamSpot) {}
0067
0068 std::unique_ptr<reco::Track> KFBasedPixelFitter::run(const std::vector<const TrackingRecHit *> &hits,
0069                                                      const TrackingRegion &region) const {
0070   std::unique_ptr<reco::Track> ret;
0071
0072   int nhits = hits.size();
0073   if (nhits < 2)
0074     return ret;
0075
0076   float ptMin = region.ptMin();
0077
0078   const GlobalPoint &vertexPos = region.origin();
0079   GlobalError vertexErr(sqr(region.originRBound()), 0, sqr(region.originRBound()), 0, 0, sqr(region.originZBound()));
0080
0081   std::vector<GlobalPoint> points(nhits);
0082   points[0] = theTracker->idToDet(hits[0]->geographicalId())->toGlobal(hits[0]->localPosition());
0083   points[1] = theTracker->idToDet(hits[1]->geographicalId())->toGlobal(hits[1]->localPosition());
0084   points[2] = theTracker->idToDet(hits[2]->geographicalId())->toGlobal(hits[2]->localPosition());
0085
0086   //
0087   //initial Kinematics
0088   //
0089   GlobalVector initMom;
0090   int charge;
0091   float theta;
0092   CircleFromThreePoints circle(points[0], points[1], points[2]);
0093   if (circle.curvature() > 1.e-4) {
0094     float invPt = PixelRecoUtilities::inversePt(circle.curvature(), *theField);
0095     float valPt = 1.f / invPt;
0096     float chargeTmp = (points[1].x() - points[0].x()) * (points[2].y() - points[1].y()) -
0097                       (points[1].y() - points[0].y()) * (points[2].x() - points[1].x());
0098     charge = (chargeTmp > 0) ? -1 : 1;
0099     float valPhi = (charge > 0) ? std::atan2(circle.center().x(), -circle.center().y())
0100                                 : std::atan2(-circle.center().x(), circle.center().y());
0101     theta = GlobalVector(points[1] - points[0]).theta();
0102     initMom = GlobalVector(valPt * cos(valPhi), valPt * sin(valPhi), valPt / tan(theta));
0103   } else {
0104     initMom = GlobalVector(points[1] - points[0]);
0105     initMom *= 10000. / initMom.perp();
0106     charge = 1;
0107     theta = initMom.theta();
0108   }
0109   GlobalTrajectoryParameters initialKine(vertexPos, initMom, TrackCharge(charge), theField);
0110
0111   //
0112   // initial error
0113   //
0114   AlgebraicSymMatrix55 C = ROOT::Math::SMatrixIdentity();
0115   float sin2th = sqr(sin(theta));
0116   float minC00 = 1.0;
0117   C[0][0] = std::max(sin2th / sqr(ptMin), minC00);
0118   float zErr = vertexErr.czz();
0119   float transverseErr = vertexErr.cxx();  // assume equal cxx cyy
0120   C[3][3] = transverseErr;
0121   C[4][4] = zErr * sin2th + transverseErr * (1 - sin2th);
0122   CurvilinearTrajectoryError initialError(C);
0123
0124   FreeTrajectoryState fts(initialKine, initialError);
0125
0126   // get updator
0127   KFUpdator updator;
0128
0129   // Now update initial state track using information from hits.
0130   TrajectoryStateOnSurface outerState;
0131   DetId outerDetId = 0;
0132   const TrackingRecHit *hit = nullptr;
0133   for (unsigned int iHit = 0; iHit < hits.size(); iHit++) {
0134     hit = hits[iHit];
0135     if (iHit == 0)
0136       outerState = thePropagator->propagate(fts, theTracker->idToDet(hit->geographicalId())->surface());
0137     outerDetId = hit->geographicalId();
0138     TrajectoryStateOnSurface state = thePropagator->propagate(outerState, theTracker->idToDet(outerDetId)->surface());
0139     if (!state.isValid())
0140       return ret;
0141     //    TransientTrackingRecHit::RecHitPointer recHit = (theTTRHBuilder->build(hit))->clone(state);
0142     TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder->build(hit);
0143     outerState = updator.update(state, *recHit);
0144     if (!outerState.isValid())
0145       return ret;
0146   }
0147
0148   TrajectoryStateOnSurface innerState = outerState;
0149   DetId innerDetId = 0;
0150   innerState.rescaleError(100000.);
0151   for (int iHit = 2; iHit >= 0; --iHit) {
0152     hit = hits[iHit];
0153     innerDetId = hit->geographicalId();
0154     TrajectoryStateOnSurface state = theOPropagator->propagate(innerState, theTracker->idToDet(innerDetId)->surface());
0155     if (!state.isValid())
0156       return ret;
0157     //  TransientTrackingRecHit::RecHitPointer recHit = (theTTRHBuilder->build(hit))->clone(state);
0158     TransientTrackingRecHit::RecHitPointer recHit = theTTRHBuilder->build(hit);
0159     innerState = updator.update(state, *recHit);
0160     if (!innerState.isValid())
0161       return ret;
0162   }
0163
0164   // extrapolate to vertex
0165   auto impactPointState = TransverseImpactPointExtrapolator(theField).extrapolate(innerState, vertexPos);
0166   if (!impactPointState.isValid())
0167     return ret;
0168
0169   //
0170   // optionally update impact point state with Bs constraint
0171   // using this potion makes sense if vertexPos (from TrackingRegion is centerewd at BeamSpot).
0172   //
0173   if (theBeamSpot) {
0174     MyBeamSpotGeomDet bsgd(Plane::build(impactPointState.surface().position(), impactPointState.surface().rotation()));
0175     MyBeamSpotHit bsrh(*theBeamSpot, &bsgd);
0176     impactPointState = updator.update(impactPointState, bsrh);  //update
0177     impactPointState =
0178         TransverseImpactPointExtrapolator(theField).extrapolate(impactPointState, vertexPos);  //reextrapolate
0179     if (!impactPointState.isValid())
0180       return ret;
0181   }
0182
0183   int ndof = 2 * hits.size() - 5;
0184   GlobalPoint vv = impactPointState.globalPosition();
0185   math::XYZPoint pos(vv.x(), vv.y(), vv.z());
0186   GlobalVector pp = impactPointState.globalMomentum();
0187   math::XYZVector mom(pp.x(), pp.y(), pp.z());
0188
0189   float chi2 = 0.;
0190   ret = std::make_unique<reco::Track>(
0191       chi2, ndof, pos, mom, impactPointState.charge(), impactPointState.curvilinearError());
0192
0193   /*
0194     vv = outerState.globalPosition();
0195     pp = outerState.globalMomentum();
0196     math::XYZPoint  outerPosition( vv.x(), vv.y(), vv.z());
0197     math::XYZVector outerMomentum( pp.x(), pp.y(), pp.z());
0198     vv = innerState.globalPosition();
0199     pp = innerState.globalMomentum();
0200     math::XYZPoint  innerPosition( vv.x(), vv.y(), vv.z());
0201     math::XYZVector innerMomentum( pp.x(), pp.y(), pp.z());
0202
0203     reco::TrackExtra extra( outerPosition, outerMomentum, true,
0204                       innerPosition, innerMomentum, true,
0205                       outerState.curvilinearError(), outerDetId,
0206                       innerState.curvilinearError(), innerDetId,
0207                       anyDirection);
0208 */
0209
0210   //  std::cout <<"TRACK CREATED" << std::endl;
0211   return ret;
0212 }