VertexFitInterface/src/TrackParameters.cc

// $Id: TrackParameters.cc,v 1.2 2008/09/27 05:48:26 loizides Exp $

#include "MitEdm/VertexFitInterface/interface/TrackParameters.h"
#include "TMath.h"
#include "MitEdm/DataFormats/interface/Types.h"
#include "MagneticField/Engine/interface/MagneticField.h"

using namespace std;
using namespace TMath;
using namespace reco;
using namespace mitedm;

//--------------------------------------------------------------------------------------------------
TrackParameters::TrackParameters(const TransientTrack &ttrk, const TrackConvention cv) :
  iConvention_(cv),
  bField_     (ttrk.field()->inTesla(GlobalPoint(0.,0.,0.)).z()),
  fCurv_      (0.5 * 0.0029979 * bField_ * -1.0)
{
  // Constructor from TransientTrack.  Track parameters are propagated back to PCA to (0,0,0)
  // in order to be consistent with the later transformation to the mvf helix
  // parameterization.  Magnetic field is also taken from the TransientTrack in this case.

  TrajectoryStateClosestToPoint traj = ttrk.trajectoryStateClosestToPoint(GlobalPoint(0.,0.,0.));
  
  GlobalPoint v = traj.theState().position();
  math::XYZPoint  pos( v.x(), v.y(), v.z() );
  GlobalVector p = traj.theState().momentum();
  math::XYZVector mom( p.x(), p.y(), p.z() );
 
  //construct intermediate reco track object, just so that we can be sure we are using the
  //correct accessors for the helix parameters
  reco::Track trk(ttrk.track().chi2(),
                  ttrk.track().ndof(),
                  pos, mom, traj.theState().charge(), 
                  traj.theState().curvilinearError(),
                  ttrk.track().algo());
  

  pars_.ResizeTo(5);
  cMat_.ResizeTo(5,5);
  if (cv == iCms) {
    pars_[0] = trk.qoverp();
    pars_[1] = trk.lambda();
    pars_[2] = trk.phi();
    pars_[3] = trk.dxy();
    pars_[4] = trk.dsz();
    for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++)
        cMat_(i,j) = trk.covariance(i,j);
    }
  }
  else {
    cout << "TrackParameters::TrackParameters -- Requested transfer: CMS -> MVF convention"
         << " not yet implemented.\n";
  }
}


//--------------------------------------------------------------------------------------------------
TrackParameters::TrackParameters(const Track* trk, const TrackConvention cv, double bField) :
  iConvention_(cv),
  bField_     (bField),
  fCurv_      (0.5 * 0.0029979 * bField_ * -1.0)
{
  // Constructor from reco::Track.  WARNING, if the helix parameters and covariance matrix stored
  // in the track are not corresponding to the PCA to (0,0,0) (NOT the case by default) then one
  // can get inconsistent results.

  pars_.ResizeTo(5);
  cMat_.ResizeTo(5,5);
  if (cv == iCms) {
    pars_[0] = trk->qoverp();
    pars_[1] = trk->lambda();
    pars_[2] = trk->phi();
    pars_[3] = trk->dxy();
    pars_[4] = trk->dsz();
    for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++)
        cMat_(i,j) = trk->covariance(i,j);
    }
  }
  else {
    cout << "TrackParameters::TrackParameters -- Requested transfer: CMS -> MVF convention"
         << " not yet implemented.\n";
  }
}

//--------------------------------------------------------------------------------------------------
TrackParameters::TrackParameters(const TrackParameters &trk) :
  iConvention_(trk.iConvention_),
  bField_     (trk.bField_),
  fCurv_      (trk.fCurv_)
{
  // Constructor.

  pars_.ResizeTo(5);
  cMat_.ResizeTo(5,5);
  for (int i=0; i<5; i++) {
    pars_(i) = trk.pars(i);
    for (int j=0; j<5; j++)
      cMat_(i,j) = trk.cMat(i,j);
  }
}

//--------------------------------------------------------------------------------------------------
TrackParameters TrackParameters::mvfTrack() const
{
  // Return a new set of track parameters in the required MVF convention.

  TrackParameters outTk;
  outTk.pars_.ResizeTo(5);
  outTk.cMat_.ResizeTo(5,5);

  if (iConvention_ == iCms) {
    outTk.setPars(0,1.0/Tan(PiOver2() - pars_[1]));  // cotTheta
    outTk.setPars(1,fCurv_*pars_[0]/Cos(pars_[1]));  // curvature
    outTk.setPars(2,pars_[4]/Cos(pars_[1]));         // z0
    outTk.setPars(3,pars_[3]);                       // d0
    double phi = pars_[2];
    while (phi<0.0)
      phi += 2*Pi();
    while (phi>2*Pi())
      phi -= 2*Pi();
    outTk.setPars(4,phi);            // phi0 [0,2pi)
    // dPidQj(i,j) gives partial dP_i/dQ_j where Pi are the new parameters and Qj are the old ones
    TMatrixD dPidQj(5,5);
    dPidQj(0,1) = 1.0/(Sin(PiOver2() - pars_[1])*Sin(PiOver2() - pars_[1]));
    dPidQj(1,0) = fCurv_/Cos(pars_[1]);
    dPidQj(1,1) = fCurv_*pars_[0]*Tan(pars_[1])/Cos(pars_[1]);
    dPidQj(2,1) = pars_[4]*Tan(pars_[1])/Cos(pars_[1]);
    dPidQj(2,4) = 1.0/Cos(pars_[1]);
    dPidQj(3,3) = 1.0;
    dPidQj(4,2) = 1.0;
    for (int i=0; i<5; i++)
      for (int j=0; j<5; j++) {
        outTk.setCMat(i,j,0.0);
        for (int k=0; k<5; k++)
          for (int l=0; l<5; l++)
            outTk.addCMat(i,j,dPidQj(i,k)*dPidQj(j,l)*cMat_(k,l));
      }
  }
  else {
    outTk.setPars(0,pars_[0]);
    outTk.setPars(1,pars_[1]);
    outTk.setPars(2,pars_[2]);
    outTk.setPars(3,pars_[3]);
    outTk.setPars(4,pars_[4]);
    for (int i=0; i<5; i++)
      for (int j=0; j<5; j++) {
        outTk.setCMat(i,j,cMat_(i,j));
      }
  }
  return outTk;
}

TrackParameters TrackParameters::cmsTrack() const
{
  // Return a new set of track parameters in the required CMS convention.

  TrackParameters outTk;
  outTk.pars_.ResizeTo(5);
  outTk.cMat_.ResizeTo(5,5);

  if (iConvention_ == iCms) {
    outTk.setPars(0,pars_[0]);
    outTk.setPars(1,pars_[1]);
    outTk.setPars(2,pars_[2]);
    outTk.setPars(3,pars_[3]);
    outTk.setPars(4,pars_[4]);
    for (int i=0; i<5; i++)
      for (int j=0; j<5; j++) {
        outTk.setCMat(i,j,cMat_(i,j));
      }
  }
  else {
    // Apply MVF to CMS conversion
    outTk.setPars(1,ATan(pars_[0]));                        // lambda
    outTk.setPars(0,pars_[1]*Cos(outTk.pars(1))/fCurv_);    // qoverp
    outTk.setPars(3,pars_[3]);                              // dxy
    outTk.setPars(4,Cos(outTk.pars(1))*pars_[2]);           // dsz
    double phi = pars_[4];
    while (phi>Pi())
      phi -= 2*Pi();                     // phi0 [-pi,pi)
    while (phi<-Pi())
      phi += 2*Pi();
    outTk.setPars(2,phi);                     // phi0 [-pi,pi)
    // dPidQj(i,j) gives partial dP_i/dQ_j where Pi are the new parameters and Qj are the old ones
    TMatrixD dPidQj(5,5);
    dPidQj(1,0) = 1.0/(1.0+pars_[0]*pars_[0]);
    dPidQj(0,0) = pars_[1]*Sin(outTk.pars(1))*dPidQj(1,0)/fCurv_;
    dPidQj(0,1) = Cos(outTk.pars(1))/fCurv_;
    dPidQj(2,4) = 1.0;
    dPidQj(3,3) = 1.0;
    dPidQj(4,0) = -pars_[2]*Sin(outTk.pars(1))*dPidQj(1,0);
    dPidQj(4,2) = Cos(outTk.pars(1));
    for (int i=0; i<5; i++)
      for (int j=0; j<5; j++) {
        outTk.setCMat(i,j,0.0);
        for (int k=0; k<5; k++)
          for (int l=0; l<5; l++)
            outTk.addCMat(i,j,dPidQj(i,k)*dPidQj(j,l)*cMat_(k,l));
      }
  }
  return outTk;
}

//--------------------------------------------------------------------------------------------------
void TrackParameters::print() const
{
  // Print track parameters.

  printf("\n==== Show track ====\n Track parameters with convention %1d.\n",int(iConvention_));
  for (int i=0; i<5; i++)
    printf(" Parameter [%1d]: %8.4f\n",i,pars_[i]);
  printf("Track covariance matrix");
  for (int i=0; i<5; i++) {
    printf("\n CV[%1d,0-4] ",i);
    for (int j=0; j<5; j++)
      printf(" %8.4f",cMat_(i,j));
  }
  printf("\n");
}
Revision:	1.3
Committed:	Mon Jan 18 14:41:34 2010 UTC (15 years, 3 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_032, Mit_031, Mit_025c_branch2, Mit_025c_branch1, Mit_030, Mit_029c, Mit_029b, Mit_030_pre1, Mit_029a, Mit_029, Mit_029_pre1, Mit_028a, Mit_025c_branch0, Mit_028, Mit_027a, Mit_027, Mit_026, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, V07-05-00, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b, ConvRejection-10-06-09, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013, Mit_013pre1, Mit_012i, HEAD
Branch point for:	Mit_025c_branch
Changes since 1.2:	+63 -8 lines
Log Message:	Fix Track Paramater phi-related bugs and add consistent propagation to origin
#	User	Rev	Content
1	bendavid	1.3	// $Id: TrackParameters.cc,v 1.2 2008/09/27 05:48:26 loizides Exp $
2	loizides	1.2
3			#include "MitEdm/VertexFitInterface/interface/TrackParameters.h"
4	paus	1.1	#include "TMath.h"
5			#include "MitEdm/DataFormats/interface/Types.h"
6	bendavid	1.3	#include "MagneticField/Engine/interface/MagneticField.h"
7	loizides	1.2
8	paus	1.1	using namespace std;
9			using namespace TMath;
10			using namespace reco;
11			using namespace mitedm;
12
13	loizides	1.2	//--------------------------------------------------------------------------------------------------
14	bendavid	1.3	TrackParameters::TrackParameters(const TransientTrack &ttrk, const TrackConvention cv) :
15			iConvention_(cv),
16			bField_ (ttrk.field()->inTesla(GlobalPoint(0.,0.,0.)).z()),
17			fCurv_ (0.5 * 0.0029979 * bField_ * -1.0)
18			{
19			// Constructor from TransientTrack. Track parameters are propagated back to PCA to (0,0,0)
20			// in order to be consistent with the later transformation to the mvf helix
21			// parameterization. Magnetic field is also taken from the TransientTrack in this case.
22
23			TrajectoryStateClosestToPoint traj = ttrk.trajectoryStateClosestToPoint(GlobalPoint(0.,0.,0.));
24
25			GlobalPoint v = traj.theState().position();
26			math::XYZPoint pos( v.x(), v.y(), v.z() );
27			GlobalVector p = traj.theState().momentum();
28			math::XYZVector mom( p.x(), p.y(), p.z() );
29
30			//construct intermediate reco track object, just so that we can be sure we are using the
31			//correct accessors for the helix parameters
32			reco::Track trk(ttrk.track().chi2(),
33			ttrk.track().ndof(),
34			pos, mom, traj.theState().charge(),
35			traj.theState().curvilinearError(),
36			ttrk.track().algo());
37
38
39			pars_.ResizeTo(5);
40			cMat_.ResizeTo(5,5);
41			if (cv == iCms) {
42			pars_[0] = trk.qoverp();
43			pars_[1] = trk.lambda();
44			pars_[2] = trk.phi();
45			pars_[3] = trk.dxy();
46			pars_[4] = trk.dsz();
47			for (int i=0; i<5; i++) {
48			for (int j=0; j<5; j++)
49			cMat_(i,j) = trk.covariance(i,j);
50			}
51			}
52			else {
53			cout << "TrackParameters::TrackParameters -- Requested transfer: CMS -> MVF convention"
54			<< " not yet implemented.\n";
55			}
56			}
57
58
59			//--------------------------------------------------------------------------------------------------
60	paus	1.1	TrackParameters::TrackParameters(const Track* trk, const TrackConvention cv, double bField) :
61			iConvention_(cv),
62			bField_ (bField),
63			fCurv_ (0.5 * 0.0029979 * bField_ * -1.0)
64			{
65	bendavid	1.3	// Constructor from reco::Track. WARNING, if the helix parameters and covariance matrix stored
66			// in the track are not corresponding to the PCA to (0,0,0) (NOT the case by default) then one
67			// can get inconsistent results.
68	loizides	1.2
69	paus	1.1	pars_.ResizeTo(5);
70			cMat_.ResizeTo(5,5);
71			if (cv == iCms) {
72			pars_[0] = trk->qoverp();
73			pars_[1] = trk->lambda();
74			pars_[2] = trk->phi();
75			pars_[3] = trk->dxy();
76			pars_[4] = trk->dsz();
77			for (int i=0; i<5; i++) {
78			for (int j=0; j<5; j++)
79			cMat_(i,j) = trk->covariance(i,j);
80			}
81			}
82			else {
83			cout << "TrackParameters::TrackParameters -- Requested transfer: CMS -> MVF convention"
84			<< " not yet implemented.\n";
85			}
86			}
87
88	loizides	1.2	//--------------------------------------------------------------------------------------------------
89	paus	1.1	TrackParameters::TrackParameters(const TrackParameters &trk) :
90			iConvention_(trk.iConvention_),
91			bField_ (trk.bField_),
92			fCurv_ (trk.fCurv_)
93			{
94	loizides	1.2	// Constructor.
95
96	paus	1.1	pars_.ResizeTo(5);
97			cMat_.ResizeTo(5,5);
98			for (int i=0; i<5; i++) {
99			pars_(i) = trk.pars(i);
100			for (int j=0; j<5; j++)
101			cMat_(i,j) = trk.cMat(i,j);
102			}
103			}
104
105	loizides	1.2	//--------------------------------------------------------------------------------------------------
106	paus	1.1	TrackParameters TrackParameters::mvfTrack() const
107			{
108	loizides	1.2	// Return a new set of track parameters in the required MVF convention.
109
110	paus	1.1	TrackParameters outTk;
111			outTk.pars_.ResizeTo(5);
112			outTk.cMat_.ResizeTo(5,5);
113
114			if (iConvention_ == iCms) {
115			outTk.setPars(0,1.0/Tan(PiOver2() - pars_[1])); // cotTheta
116			outTk.setPars(1,fCurv_*pars_[0]/Cos(pars_[1])); // curvature
117			outTk.setPars(2,pars_[4]/Cos(pars_[1])); // z0
118			outTk.setPars(3,pars_[3]); // d0
119	bendavid	1.3	double phi = pars_[2];
120			while (phi<0.0)
121			phi += 2*Pi();
122			while (phi>2*Pi())
123			phi -= 2*Pi();
124			outTk.setPars(4,phi); // phi0 [0,2pi)
125	paus	1.1	// dPidQj(i,j) gives partial dP_i/dQ_j where Pi are the new parameters and Qj are the old ones
126			TMatrixD dPidQj(5,5);
127			dPidQj(0,1) = 1.0/(Sin(PiOver2() - pars_[1])*Sin(PiOver2() - pars_[1]));
128			dPidQj(1,0) = fCurv_/Cos(pars_[1]);
129			dPidQj(1,1) = fCurv_pars_[0]Tan(pars_[1])/Cos(pars_[1]);
130			dPidQj(2,1) = pars_[4]*Tan(pars_[1])/Cos(pars_[1]);
131			dPidQj(2,4) = 1.0/Cos(pars_[1]);
132			dPidQj(3,3) = 1.0;
133			dPidQj(4,2) = 1.0;
134			for (int i=0; i<5; i++)
135			for (int j=0; j<5; j++) {
136			outTk.setCMat(i,j,0.0);
137			for (int k=0; k<5; k++)
138			for (int l=0; l<5; l++)
139			outTk.addCMat(i,j,dPidQj(i,k)dPidQj(j,l)cMat_(k,l));
140			}
141			}
142			else {
143			outTk.setPars(0,pars_[0]);
144			outTk.setPars(1,pars_[1]);
145			outTk.setPars(2,pars_[2]);
146			outTk.setPars(3,pars_[3]);
147			outTk.setPars(4,pars_[4]);
148			for (int i=0; i<5; i++)
149			for (int j=0; j<5; j++) {
150			outTk.setCMat(i,j,cMat_(i,j));
151			}
152			}
153			return outTk;
154			}
155
156			TrackParameters TrackParameters::cmsTrack() const
157			{
158	loizides	1.2	// Return a new set of track parameters in the required CMS convention.
159
160	paus	1.1	TrackParameters outTk;
161			outTk.pars_.ResizeTo(5);
162			outTk.cMat_.ResizeTo(5,5);
163
164			if (iConvention_ == iCms) {
165			outTk.setPars(0,pars_[0]);
166			outTk.setPars(1,pars_[1]);
167			outTk.setPars(2,pars_[2]);
168			outTk.setPars(3,pars_[3]);
169			outTk.setPars(4,pars_[4]);
170			for (int i=0; i<5; i++)
171			for (int j=0; j<5; j++) {
172			outTk.setCMat(i,j,cMat_(i,j));
173			}
174			}
175			else {
176			// Apply MVF to CMS conversion
177			outTk.setPars(1,ATan(pars_[0])); // lambda
178			outTk.setPars(0,pars_[1]*Cos(outTk.pars(1))/fCurv_); // qoverp
179			outTk.setPars(3,pars_[3]); // dxy
180			outTk.setPars(4,Cos(outTk.pars(1))*pars_[2]); // dsz
181	bendavid	1.3	double phi = pars_[4];
182			while (phi>Pi())
183			phi -= 2*Pi(); // phi0 [-pi,pi)
184			while (phi<-Pi())
185			phi += 2*Pi();
186			outTk.setPars(2,phi); // phi0 [-pi,pi)
187	paus	1.1	// dPidQj(i,j) gives partial dP_i/dQ_j where Pi are the new parameters and Qj are the old ones
188			TMatrixD dPidQj(5,5);
189			dPidQj(1,0) = 1.0/(1.0+pars_[0]*pars_[0]);
190			dPidQj(0,0) = pars_[1]Sin(outTk.pars(1))dPidQj(1,0)/fCurv_;
191			dPidQj(0,1) = Cos(outTk.pars(1))/fCurv_;
192			dPidQj(2,4) = 1.0;
193			dPidQj(3,3) = 1.0;
194			dPidQj(4,0) = -pars_[2]Sin(outTk.pars(1))dPidQj(1,0);
195			dPidQj(4,2) = Cos(outTk.pars(1));
196			for (int i=0; i<5; i++)
197			for (int j=0; j<5; j++) {
198			outTk.setCMat(i,j,0.0);
199			for (int k=0; k<5; k++)
200			for (int l=0; l<5; l++)
201			outTk.addCMat(i,j,dPidQj(i,k)dPidQj(j,l)cMat_(k,l));
202			}
203			}
204			return outTk;
205			}
206
207	loizides	1.2	//--------------------------------------------------------------------------------------------------
208	paus	1.1	void TrackParameters::print() const
209			{
210	loizides	1.2	// Print track parameters.
211
212	paus	1.1	printf("\n==== Show track ====\n Track parameters with convention %1d.\n",int(iConvention_));
213			for (int i=0; i<5; i++)
214			printf(" Parameter [%1d]: %8.4f\n",i,pars_[i]);
215			printf("Track covariance matrix");
216			for (int i=0; i<5; i++) {
217			printf("\n CV[%1d,0-4] ",i);
218			for (int j=0; j<5; j++)
219			printf(" %8.4f",cMat_(i,j));
220			}
221			printf("\n");
222			}