PixelTrackFitting/src/PixelNtupletsFitter.cc

0001 #include "CommonTools/Utils/interface/DynArray.h"
0002 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0003 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0004 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0005 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0006 #include "DataFormats/GeometryVector/interface/Pi.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/GeomDetType.h"
0012 #include "MagneticField/Engine/interface/MagneticField.h"
0013 #include "RecoTracker/PixelTrackFitting/interface/BrokenLine.h"
0014 #include "RecoTracker/PixelTrackFitting/interface/PixelNtupletsFitter.h"
0015 #include "RecoTracker/PixelTrackFitting/interface/PixelTrackBuilder.h"
0016 #include "RecoTracker/PixelTrackFitting/interface/PixelTrackErrorParam.h"
0017 #include "RecoTracker/PixelTrackFitting/interface/RiemannFit.h"
0018 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0019
0020 using namespace std;
0021
0022 PixelNtupletsFitter::PixelNtupletsFitter(float nominalB, const MagneticField* field, bool useRiemannFit)
0023     : nominalB_(nominalB), field_(field), useRiemannFit_(useRiemannFit) {}
0024
0025 std::unique_ptr<reco::Track> PixelNtupletsFitter::run(const std::vector<const TrackingRecHit*>& hits,
0026                                                       const TrackingRegion& region) const {
0027   using namespace riemannFit;
0028
0029   std::unique_ptr<reco::Track> ret;
0030
0031   unsigned int nhits = hits.size();
0032
0033   if (nhits < 2)
0034     return ret;
0035
0036   declareDynArray(GlobalPoint, nhits, points);
0037   declareDynArray(GlobalError, nhits, errors);
0038   declareDynArray(bool, nhits, isBarrel);
0039
0040   for (unsigned int i = 0; i != nhits; ++i) {
0041     auto const& recHit = hits[i];
0042     points[i] = GlobalPoint(recHit->globalPosition().basicVector() - region.origin().basicVector());
0043     errors[i] = recHit->globalPositionError();
0044     isBarrel[i] = recHit->detUnit()->type().isBarrel();
0045   }
0046
0047   assert(nhits == 4);
0048   riemannFit::Matrix3xNd<4> hits_gp;
0049
0050   Eigen::Matrix<float, 6, 4> hits_ge = Eigen::Matrix<float, 6, 4>::Zero();
0051
0052   for (unsigned int i = 0; i < nhits; ++i) {
0053     hits_gp.col(i) << points[i].x(), points[i].y(), points[i].z();
0054
0055     hits_ge.col(i) << errors[i].cxx(), errors[i].cyx(), errors[i].cyy(), errors[i].czx(), errors[i].czy(),
0056         errors[i].czz();
0057   }
0058
0059   HelixFit fittedTrack = useRiemannFit_ ? riemannFit::helixFit(hits_gp, hits_ge, nominalB_, true)
0060                                         : brokenline::helixFit(hits_gp, hits_ge, nominalB_);
0061
0062   int iCharge = fittedTrack.qCharge;
0063
0064   // parameters are:
0065   // 0: phi
0066   // 1: tip
0067   // 2: curvature
0068   // 3: cottheta
0069   // 4: zip
0070   float valPhi = fittedTrack.par(0);
0071
0072   float valTip = fittedTrack.par(1);
0073
0074   float valCotTheta = fittedTrack.par(3);
0075
0076   float valZip = fittedTrack.par(4);
0077   float valPt = fittedTrack.par(2);
0078   //
0079   //  PixelTrackErrorParam param(valEta, valPt);
0080   float errValPhi = std::sqrt(fittedTrack.cov(0, 0));
0081   float errValTip = std::sqrt(fittedTrack.cov(1, 1));
0082
0083   float errValPt = std::sqrt(fittedTrack.cov(2, 2));
0084
0085   float errValCotTheta = std::sqrt(fittedTrack.cov(3, 3));
0086   float errValZip = std::sqrt(fittedTrack.cov(4, 4));
0087
0088   float chi2 = fittedTrack.chi2_line + fittedTrack.chi2_circle;
0089
0090   PixelTrackBuilder builder;
0091   Measurement1D phi(valPhi, errValPhi);
0092   Measurement1D tip(valTip, errValTip);
0093
0094   Measurement1D pt(valPt, errValPt);
0095   Measurement1D cotTheta(valCotTheta, errValCotTheta);
0096   Measurement1D zip(valZip, errValZip);
0097
0098   ret.reset(builder.build(pt, phi, cotTheta, tip, zip, chi2, iCharge, hits, field_, region.origin()));
0099   return ret;
0100 }