CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoPixelVertexing/PixelTrackFitting/src/PixelTrackBuilder.cc

Go to the documentation of this file.
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     // TrajectoryStateOnSurface state(bstate);
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   // BTSOS hold BP in a shared pointer and  will be autodeleted when BTSOS goes out of scope...
00179   // to avoid memory churn we allocate it locally and just avoid it be deleted by refcount... 
00180   Plane impPointPlane(origin, rot);
00181   // (twice just to be sure!)
00182   impPointPlane.addReference(); impPointPlane.addReference();
00183   // use Base (too avoid a useless new) 
00184   BasicTrajectoryStateOnSurface impactPointState( lpar , error, impPointPlane, mf, 1.0);
00185   
00186   //checkState(impactPointState,mf);
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