00001 #include "PixelTrackBuilder.h"
00002
00003
00004 #include "DataFormats/GeometrySurface/interface/LocalError.h"
00005 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
00006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00007
00008 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryError.h"
00009 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
00010 #include "TrackingTools/TrajectoryState/interface/BasicTrajectoryStateOnSurface.h"
00011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00012
00013 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
00014 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
00015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00016 #include "DataFormats/TrackReco/interface/TrackBase.h"
00017
00018 #include <sstream>
00019 using namespace std;
00020 using namespace reco;
00021
00022
00023 template <class T> T inline sqr( T t) {return t*t;}
00024
00025 namespace {
00026 std::string print(
00027 const Measurement1D & pt,
00028 const Measurement1D & phi,
00029 const Measurement1D & cotTheta,
00030 const Measurement1D & tip,
00031 const Measurement1D & zip,
00032 float chi2,
00033 int charge)
00034 {
00035 ostringstream str;
00036 str <<"\t pt: " << pt.value() <<"+/-"<<pt.error()
00037 <<"\t phi: " << phi.value() <<"+/-"<<phi.error()
00038 <<"\t cot: " << cotTheta.value() <<"+/-"<<cotTheta.error()
00039 <<"\t tip: " << tip.value() <<"+/-"<<tip.error()
00040 <<"\t zip: " << zip.value() <<"+/-"<<zip.error()
00041 <<"\t charge: " << charge;
00042 return str.str();
00043 }
00044
00045 std::string print(const reco::Track & track, const GlobalPoint & origin)
00046 {
00047
00048 math::XYZPoint bs(origin.x(), origin.y(), origin.z());
00049
00050 Measurement1D phi( track.phi(), track.phiError());
00051
00052 float theta = track.theta();
00053 float cosTheta = cos(theta);
00054 float sinTheta = sin(theta);
00055 float errLambda2 = sqr( track.lambdaError() );
00056 Measurement1D cotTheta(cosTheta/sinTheta, sqrt(errLambda2)/sqr(sinTheta));
00057
00058 float pt_v = track.pt();
00059 float errInvP2 = sqr(track.qoverpError());
00060 float covIPtTheta = track.covariance(TrackBase::i_qoverp, TrackBase::i_lambda);
00061 float errInvPt2 = ( errInvP2
00062 + sqr(cosTheta/pt_v)*errLambda2
00063 + 2*(cosTheta/pt_v)*covIPtTheta
00064 ) / sqr(sinTheta);
00065 Measurement1D pt(pt_v, sqr(pt_v)*sqrt(errInvPt2));
00066
00067 Measurement1D tip(track.dxy(bs), track.d0Error());
00068
00069 Measurement1D zip(track.dz(bs), track.dzError());
00070
00071 return print(pt, phi, cotTheta, tip, zip, track.chi2(), track.charge());
00072 }
00073
00074 std::string print(const BasicTrajectoryStateOnSurface & state)
00075 {
00076
00077 float pt_v = state.globalMomentum().perp();
00078 float phi_v = state.globalMomentum().phi();
00079 float theta_v = state.globalMomentum().theta();
00080
00081 CurvilinearTrajectoryError curv = state.curvilinearError();
00082 float errPhi2 = curv.matrix()(3,3);
00083 float errLambda2 = curv.matrix()(2,2);
00084 float errInvP2 = curv.matrix()(1,1);
00085 float covIPtTheta = curv.matrix()(1,2);
00086 float cosTheta = cos(theta_v);
00087 float sinTheta = sin(theta_v);
00088 float errInvPt2 = ( errInvP2
00089 + sqr(cosTheta/pt_v)*errLambda2
00090 + 2*(cosTheta/pt_v)*covIPtTheta) / sqr(sinTheta);
00091 float errCotTheta = sqrt(errLambda2)/sqr(sinTheta) ;
00092 Measurement1D pt(pt_v, sqr(pt_v) * sqrt(errInvPt2));
00093 Measurement1D phi(phi_v, sqrt(errPhi2) );
00094 Measurement1D cotTheta(cosTheta/sinTheta, errCotTheta);
00095
00096 float zip_v = state.globalPosition().z();
00097 float zip_e = sqrt( state.localError().matrix()(4,4));
00098 Measurement1D zip(zip_v, zip_e);
00099
00100 float tip_v = state.localPosition().x();
00101 int tip_sign = (state.localMomentum().y()*cotTheta.value() > 0) ? -1 : 1;
00102 float tip_e = sqrt( state.localError().matrix()(3,3) );
00103 Measurement1D tip( tip_sign*tip_v, tip_e);
00104
00105 return print(pt, phi, cotTheta, tip, zip, 0., state.charge());
00106 }
00107
00108
00109 void checkState(const BasicTrajectoryStateOnSurface & bstate, const MagneticField* mf, const GlobalPoint & origin)
00110 {
00111 TrajectoryStateOnSurface state(bstate.clone());
00112
00113 LogTrace("")<<" *** PixelTrackBuilder::checkState: ";
00114 LogTrace("")<<"INPUT, ROTATION" << endl<<state.surface().rotation();
00115 LogTrace("")<<"INPUT, TSOS:"<<endl<<state;
00116
00117 TransverseImpactPointExtrapolator tipe(mf);
00118 TrajectoryStateOnSurface test= tipe.extrapolate(state, origin);
00119 LogTrace("")<<"CHECK-1 ROTATION" << endl<<"\n"<<test.surface().rotation();
00120 LogTrace("")<<"CHECK-1 TSOS" << endl<<test;
00121
00122 TSCPBuilderNoMaterial tscpBuilder;
00123 TrajectoryStateClosestToPoint tscp =
00124 tscpBuilder(*(state.freeState()), origin);
00125 FreeTrajectoryState fs = tscp.theState();
00126 LogTrace("")<<"CHECK-2 FTS: " << fs;
00127 }
00128
00129 }
00130
00131
00132 reco::Track * PixelTrackBuilder::build(
00133 const Measurement1D & pt,
00134 const Measurement1D & phi,
00135 const Measurement1D & cotTheta,
00136 const Measurement1D & tip,
00137 const Measurement1D & zip,
00138 float chi2,
00139 int charge,
00140 const std::vector<const TrackingRecHit* >& hits,
00141 const MagneticField * mf,
00142 const GlobalPoint & origin) const
00143 {
00144
00145 LogDebug("PixelTrackBuilder::build");
00146 LogTrace("")<<"reconstructed TRIPLET kinematics:\n"<<print(pt,phi,cotTheta,tip,zip,chi2,charge);
00147
00148 double sinTheta = 1/std::sqrt(1+sqr(cotTheta.value()));
00149 double cosTheta = cotTheta.value()*sinTheta;
00150 int tipSign = tip.value() > 0 ? 1 : -1;
00151
00152 AlgebraicSymMatrix55 m;
00153 double invPtErr = 1./sqr(pt.value()) * pt.error();
00154 m(0,0) = sqr(sinTheta) * (
00155 sqr(invPtErr)
00156 + sqr(cotTheta.error()/pt.value()*cosTheta * sinTheta)
00157 );
00158 m(0,2) = sqr( cotTheta.error()) * cosTheta * sqr(sinTheta) / pt.value();
00159 m(1,1) = sqr( phi.error() );
00160 m(2,2) = sqr( cotTheta.error());
00161 m(3,3) = sqr( tip.error() );
00162 m(4,4) = sqr( zip.error() );
00163 LocalTrajectoryError error(m);
00164
00165 LocalTrajectoryParameters lpar(
00166 LocalPoint(tipSign*tip.value(), -tipSign*zip.value(), 0),
00167 LocalVector(0., -tipSign*pt.value()*cotTheta.value(), pt.value()),
00168 charge);
00169
00170
00171 float sp = std::sin(phi.value());
00172 float cp = std::cos(phi.value());
00173 Surface::RotationType rot(
00174 sp*tipSign, -cp*tipSign, 0,
00175 0 , 0, -tipSign,
00176 cp , sp , 0);
00177
00178
00179
00180 Plane impPointPlane(origin, rot);
00181
00182 impPointPlane.addReference(); impPointPlane.addReference();
00183
00184 BasicTrajectoryStateOnSurface impactPointState( lpar , error, impPointPlane, mf, 1.0);
00185
00186
00187 LogTrace("")<<"constructed TSOS :\n"<<print(impactPointState);
00188
00189 int ndof = 2*hits.size()-5;
00190 GlobalPoint vv = impactPointState.globalPosition();
00191 math::XYZPoint pos( vv.x(), vv.y(), vv.z() );
00192 GlobalVector pp = impactPointState.globalMomentum();
00193 math::XYZVector mom( pp.x(), pp.y(), pp.z() );
00194
00195 reco::Track * track = new reco::Track( chi2, ndof, pos, mom,
00196 impactPointState.charge(), impactPointState.curvilinearError());
00197
00198 LogTrace("") <<"RECONSTRUCTED TRACK (0,0,0):\n"<< print(*track,GlobalPoint(0,0,0))<<std::endl;;
00199 LogTrace("") <<"RECONSTRUCTED TRACK "<<origin<<"\n"<< print(*track,origin)<<std::endl;;
00200
00201 return track;
00202 }
00203