PixelTrackFitting/src/PixelTrackBuilder.cc

0001 #include "RecoTracker/PixelTrackFitting/interface/PixelTrackBuilder.h"
0002
0003 #include "DataFormats/GeometrySurface/interface/LocalError.h"
0004 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0005 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0006
0007 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryError.h"
0008 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
0009 #include "TrackingTools/TrajectoryState/interface/BasicTrajectoryStateOnSurface.h"
0010 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0011
0012 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0013 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "DataFormats/TrackReco/interface/TrackBase.h"
0016
0017 #include <sstream>
0018 using namespace std;
0019 using namespace reco;
0020
0021 template <class T>
0022 T inline sqr(T t) {
0023   return t * t;
0024 }
0025
0026 namespace {
0027   __attribute__((unused)) std::string print(const Measurement1D& pt,
0028                                             const Measurement1D& phi,
0029                                             const Measurement1D& cotTheta,
0030                                             const Measurement1D& tip,
0031                                             const Measurement1D& zip,
0032                                             float chi2,
0033                                             int charge) {
0034     ostringstream str;
0035     str << "\t pt: " << pt.value() << "+/-" << pt.error() << "\t phi: " << phi.value() << "+/-" << phi.error()
0036         << "\t cot: " << cotTheta.value() << "+/-" << cotTheta.error() << "\t tip: " << tip.value() << "+/-"
0037         << tip.error() << "\t zip: " << zip.value() << "+/-" << zip.error() << "\t chi2: " << chi2
0038         << "\t charge: " << charge;
0039     return str.str();
0040   }
0041
0042   __attribute__((unused)) std::string print(const reco::Track& track, const GlobalPoint& origin) {
0043     math::XYZPoint bs(origin.x(), origin.y(), origin.z());
0044
0045     Measurement1D phi(track.phi(), track.phiError());
0046
0047     float theta = track.theta();
0048     float cosTheta = cos(theta);
0049     float sinTheta = sin(theta);
0050     float errLambda2 = sqr(track.lambdaError());
0051     Measurement1D cotTheta(cosTheta / sinTheta, sqrt(errLambda2) / sqr(sinTheta));
0052
0053     float pt_v = track.pt();
0054     float errInvP2 = sqr(track.qoverpError());
0055     float covIPtTheta = track.covariance(TrackBase::i_qoverp, TrackBase::i_lambda);
0056     float errInvPt2 =
0057         (errInvP2 + sqr(cosTheta / pt_v) * errLambda2 + 2 * (cosTheta / pt_v) * covIPtTheta) / sqr(sinTheta);
0058     Measurement1D pt(pt_v, sqr(pt_v) * sqrt(errInvPt2));
0059
0060     Measurement1D tip(track.dxy(bs), track.d0Error());
0061
0062     Measurement1D zip(track.dz(bs), track.dzError());
0063
0064     return print(pt, phi, cotTheta, tip, zip, track.chi2(), track.charge());
0065   }
0066
0067   __attribute__((unused)) std::string print(const BasicTrajectoryStateOnSurface& state) {
0068     // TrajectoryStateOnSurface state(bstate);
0069     float pt_v = state.globalMomentum().perp();
0070     float phi_v = state.globalMomentum().phi();
0071     float theta_v = state.globalMomentum().theta();
0072
0073     CurvilinearTrajectoryError curv = state.curvilinearError();
0074     float errPhi2 = curv.matrix()(3, 3);
0075     float errLambda2 = curv.matrix()(2, 2);
0076     float errInvP2 = curv.matrix()(1, 1);
0077     float covIPtTheta = curv.matrix()(1, 2);
0078     float cosTheta = cos(theta_v);
0079     float sinTheta = sin(theta_v);
0080     float errInvPt2 =
0081         (errInvP2 + sqr(cosTheta / pt_v) * errLambda2 + 2 * (cosTheta / pt_v) * covIPtTheta) / sqr(sinTheta);
0082     float errCotTheta = sqrt(errLambda2) / sqr(sinTheta);
0083     Measurement1D pt(pt_v, sqr(pt_v) * sqrt(errInvPt2));
0084     Measurement1D phi(phi_v, sqrt(errPhi2));
0085     Measurement1D cotTheta(cosTheta / sinTheta, errCotTheta);
0086
0087     float zip_v = state.globalPosition().z();
0088     float zip_e = sqrt(state.localError().matrix()(4, 4));
0089     Measurement1D zip(zip_v, zip_e);
0090
0091     float tip_v = state.localPosition().x();
0092     int tip_sign = (state.localMomentum().y() * cotTheta.value() > 0) ? -1 : 1;
0093     float tip_e = sqrt(state.localError().matrix()(3, 3));
0094     Measurement1D tip(tip_sign * tip_v, tip_e);
0095
0096     return print(pt, phi, cotTheta, tip, zip, 0., state.charge());
0097   }
0098
0099 #ifdef DEBUG_STATE
0100   inline void checkState(const BasicTrajectoryStateOnSurface& bstate,
0101                          const MagneticField* mf,
0102                          const GlobalPoint& origin) {
0103     TrajectoryStateOnSurface state(bstate.clone());
0104
0105     LogTrace("") << " *** PixelTrackBuilder::checkState: ";
0106     LogTrace("") << "INPUT,  ROTATION" << '\n' << state.surface().rotation();
0107     LogTrace("") << "INPUT,  TSOS:" << '\n' << state;
0108
0109     TransverseImpactPointExtrapolator tipe(mf);
0110     TrajectoryStateOnSurface test = tipe.extrapolate(state, origin);
0111     LogTrace("") << "CHECK-1 ROTATION" << '\n' << test.surface().rotation();
0112     LogTrace("") << "CHECK-1 TSOS" << '\n' << test;
0113
0114     TSCPBuilderNoMaterial tscpBuilder;
0115     TrajectoryStateClosestToPoint tscp = tscpBuilder(*(state.freeState()), origin);
0116     const FreeTrajectoryState& fs = tscp.theState();
0117     LogTrace("") << "CHECK-2 FTS: " << fs;
0118   }
0119 #endif
0120
0121 }  // namespace
0122
0123 reco::Track* PixelTrackBuilder::build(const Measurement1D& pt,
0124                                       const Measurement1D& phi,
0125                                       const Measurement1D& cotTheta,
0126                                       const Measurement1D& tip,
0127                                       const Measurement1D& zip,
0128                                       float chi2,
0129                                       int charge,
0130                                       const std::vector<const TrackingRecHit*>& hits,
0131                                       const MagneticField* mf,
0132                                       const GlobalPoint& origin) const {
0133   LogDebug("PixelTrackBuilder::build") << "";
0134   LogTrace("") << "Reconstructed triplet kinematics: " << print(pt, phi, cotTheta, tip, zip, chi2, charge);
0135
0136   double sinTheta = 1 / std::sqrt(1 + sqr(cotTheta.value()));
0137   double cosTheta = cotTheta.value() * sinTheta;
0138   int tipSign = tip.value() > 0 ? 1 : -1;
0139
0140   AlgebraicSymMatrix55 m;
0141   double invPtErr = 1. / sqr(pt.value()) * pt.error();
0142   m(0, 0) = sqr(sinTheta) * (sqr(invPtErr) + sqr(cotTheta.error() / pt.value() * cosTheta * sinTheta));
0143   m(0, 2) = sqr(cotTheta.error()) * cosTheta * sqr(sinTheta) / pt.value();
0144   m(1, 1) = sqr(phi.error());
0145   m(2, 2) = sqr(cotTheta.error());
0146   m(3, 3) = sqr(tip.error());
0147   m(4, 4) = sqr(zip.error());
0148   LocalTrajectoryError error(m);
0149
0150   LocalTrajectoryParameters lpar(LocalPoint(tipSign * tip.value(), -tipSign * zip.value(), 0),
0151                                  LocalVector(0., -tipSign * pt.value() * cotTheta.value(), pt.value()),
0152                                  charge);
0153
0154   float sp = std::sin(phi.value());
0155   float cp = std::cos(phi.value());
0156   Surface::RotationType rot(sp * tipSign, -cp * tipSign, 0, 0, 0, -tipSign, cp, sp, 0);
0157
0158   // BTSOS hold Surface in a shared pointer and  will be autodeleted when BTSOS goes out of scope...
0159   // to avoid memory churn we allocate it locally and just avoid it be deleted by refcount...
0160   Plane impPointPlane(origin, rot);
0161   // (twice just to be sure!)
0162   impPointPlane.addReference();
0163   impPointPlane.addReference();
0164   // use Base (to avoid a useless new)
0165   BasicTrajectoryStateOnSurface impactPointState(lpar, error, impPointPlane, mf);
0166
0167 #ifdef DEBUG_STATE
0168   checkState(impactPointState, mf);
0169 #endif
0170   LogTrace("") << "constructed TSOS:\n" << print(impactPointState);
0171
0172   int ndof = 2 * hits.size() - 5;
0173   GlobalPoint vv = impactPointState.globalPosition();
0174   math::XYZPoint pos(vv.x(), vv.y(), vv.z());
0175   GlobalVector pp = impactPointState.globalMomentum();
0176   math::XYZVector mom(pp.x(), pp.y(), pp.z());
0177
0178   reco::Track* track =
0179       new reco::Track(chi2, ndof, pos, mom, impactPointState.charge(), impactPointState.curvilinearError());
0180
0181   return track;
0182 }