CMSSW Reference Manual

AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian	(	const GlobalTrajectoryParameters &	globalParameters,
		const GlobalPoint &	x,
		const GlobalVector &	p,
		const double &	s
	)

get Field at starting state (internally)

Definition at line 5 of file AnalyticalCurvilinearJacobian.cc.

References e, h, GlobalTrajectoryParameters::magneticFieldInInverseGeV(), GlobalTrajectoryParameters::position(), and GlobalTrajectoryParameters::transverseCurvature().

00008                   : theJacobian(AlgebraicMatrixID()) 
00009 {
00010   //
00011   // helix: calculate full jacobian
00012   //
00013   if ( s*s*fabs(globalParameters.transverseCurvature())>1.e-5 ) { 
00014     GlobalPoint xStart = globalParameters.position();
00015     GlobalVector h  = globalParameters.magneticFieldInInverseGeV(xStart);
00016     computeFullJacobian(globalParameters,x,p,h,s);
00017   }
00018   //
00019   // straight line approximation, error in RPhi about 0.1um
00020   //
00021   else
00022     computeStraightLineJacobian(globalParameters,x,p,s);
00023    //dbg::dbg_trace(1,"ACJ1", globalParameters.vector(),x,p,s,theJacobian);
00024 }

AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian	(	const GlobalTrajectoryParameters &	globalParameters,
		const GlobalPoint &	x,
		const GlobalVector &	p,
		const GlobalVector &	theFieldInInverseGeV,
		const double &	s
	)

new: give Field as a parameter

Definition at line 28 of file AnalyticalCurvilinearJacobian.cc.

References e, and GlobalTrajectoryParameters::transverseCurvature().

00032                   : theJacobian(AlgebraicMatrixID()) 
00033 {
00034   //
00035   // helix: calculate full jacobian
00036   //
00037   if ( s*s*fabs(globalParameters.transverseCurvature())>1.e-5 )
00038     computeFullJacobian(globalParameters,x,p,h,s);
00039   //
00040   // straight line approximation, error in RPhi about 0.1um
00041   //
00042   else
00043     computeStraightLineJacobian(globalParameters,x,p,s);
00044   
00045   //dbg::dbg_trace(1,"ACJ2", globalParameters.vector(),x,p,s,theJacobian);
00046 }

virtual AnalyticalCurvilinearJacobian::~AnalyticalCurvilinearJacobian ( ) [inline, virtual]

Definition at line 36 of file AnalyticalCurvilinearJacobian.h.

00036 {}

void AnalyticalCurvilinearJacobian::computeFullJacobian	(	const GlobalTrajectoryParameters &	globalParameters,
		const GlobalPoint &	x,
		const GlobalVector &	p,
		const GlobalVector &	h,
		const double &	s
	)			`[private]`

result for non-vanishing curvature

Definition at line 51 of file AnalyticalCurvilinearJacobian.cc.

References funct::cos(), prof2calltree::cost, h1, h2, h3, PV3DBase< T, PVType, FrameType >::mag(), GlobalTrajectoryParameters::momentum(), p1, p2, PV3DBase< T, PVType, FrameType >::perp(), GlobalTrajectoryParameters::position(), s2, s3, GlobalTrajectoryParameters::signedInverseMomentum(), funct::sin(), funct::sqrt(), theta, Vector3DBase< T, FrameTag >::unit(), v11, v12, v13, v21, v22, v23, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

00056 {    
00057   //GlobalVector p1 = fts.momentum().unit();
00058   GlobalVector p1 = globalParameters.momentum().unit();
00059   GlobalVector p2 = p.unit();
00060   //GlobalPoint xStart = fts.position();
00061   GlobalPoint xStart = globalParameters.position();
00062   GlobalVector dx = xStart - x;
00063   //GlobalVector h  = MagneticField::inInverseGeV(xStart);
00064   // Martijn: field is now given as parameter.. GlobalVector h  = globalParameters.magneticFieldInInverseGeV(xStart);
00065 
00066   //double qbp = fts.signedInverseMomentum();
00067   double qbp = globalParameters.signedInverseMomentum();
00068   double absS = s;
00069   
00070   // calculate transport matrix
00071   // Origin: TRPRFN
00072   double t11 = p1.x(); double t12 = p1.y(); double t13 = p1.z();
00073   double t21 = p2.x(); double t22 = p2.y(); double t23 = p2.z();
00074   double cosl0 = p1.perp(); double cosl1 = 1./p2.perp();
00075   //AlgebraicMatrix a(5,5,1);
00076   // define average magnetic field and gradient 
00077   // at initial point - inlike TRPRFN
00078   GlobalVector hn = h.unit();
00079   double qp = -h.mag();
00080 //   double q = -h.mag()*qbp;
00081   double q = qp*qbp;
00082   double theta = q*absS; double sint = sin(theta); double cost = cos(theta);
00083   double hn1 = hn.x(); double hn2 = hn.y(); double hn3 = hn.z();
00084   double dx1 = dx.x(); double dx2 = dx.y(); double dx3 = dx.z();
00085   double gamma = hn1*t21 + hn2*t22 + hn3*t23;
00086   double an1 = hn2*t23 - hn3*t22;
00087   double an2 = hn3*t21 - hn1*t23;
00088   double an3 = hn1*t22 - hn2*t21;
00089   double au = 1./sqrt(t11*t11 + t12*t12);
00090   double u11 = -au*t12; double u12 = au*t11;
00091   double v11 = -t13*u12; double v12 = t13*u11; double v13 = t11*u12 - t12*u11;
00092   au = 1./sqrt(t21*t21 + t22*t22);
00093   double u21 = -au*t22; double u22 = au*t21;
00094   double v21 = -t23*u22; double v22 = t23*u21; double v23 = t21*u22 - t22*u21;
00095   // now prepare the transport matrix
00096   // pp only needed in high-p case (WA)
00097 //   double pp = 1./qbp;
00099 // moved up (where -h.mag() is needed()
00100 //   double qp = q*pp;
00101   double anv = -(hn1*u21 + hn2*u22          );
00102   double anu =  (hn1*v21 + hn2*v22 + hn3*v23);
00103   double omcost = 1. - cost; double tmsint = theta - sint;
00104   
00105   double hu1 =         - hn3*u12;
00106   double hu2 = hn3*u11;
00107   double hu3 = hn1*u12 - hn2*u11;
00108   
00109   double hv1 = hn2*v13 - hn3*v12;
00110   double hv2 = hn3*v11 - hn1*v13;
00111   double hv3 = hn1*v12 - hn2*v11;
00112   
00113   //   1/p - doesn't change since |p1| = |p2|
00114   
00115   //   lambda
00116   
00117   theJacobian(1,0) = -qp*anv*(t21*dx1 + t22*dx2 + t23*dx3);
00118   
00119   theJacobian(1,1) = cost*(v11*v21 + v12*v22 + v13*v23) +
00120                       sint*(hv1*v21 + hv2*v22 + hv3*v23) +
00121                     omcost*(hn1*v11 + hn2*v12 + hn3*v13) *
00122                            (hn1*v21 + hn2*v22 + hn3*v23) +
00123                 anv*(-sint*(v11*t21 + v12*t22 + v13*t23) +
00124                     omcost*(v11*an1 + v12*an2 + v13*an3) -
00125               tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
00126 
00127   theJacobian(1,2) = cost*(u11*v21 + u12*v22          ) +
00128                       sint*(hu1*v21 + hu2*v22 + hu3*v23) +
00129                     omcost*(hn1*u11 + hn2*u12          ) *
00130                            (hn1*v21 + hn2*v22 + hn3*v23) +
00131                 anv*(-sint*(u11*t21 + u12*t22          ) +
00132                     omcost*(u11*an1 + u12*an2          ) -
00133               tmsint*gamma*(hn1*u11 + hn2*u12          ) );
00134   theJacobian(1,2) *= cosl0;
00135 
00136   theJacobian(1,3) = -q*anv*(u11*t21 + u12*t22          );
00137 
00138   theJacobian(1,4) = -q*anv*(v11*t21 + v12*t22 + v13*t23);
00139 
00140   //   phi
00141 
00142   theJacobian(2,0) = -qp*anu*(t21*dx1 + t22*dx2 + t23*dx3)*cosl1;
00143 
00144   theJacobian(2,1) = cost*(v11*u21 + v12*u22          ) +
00145                       sint*(hv1*u21 + hv2*u22          ) +
00146                     omcost*(hn1*v11 + hn2*v12 + hn3*v13) *
00147                            (hn1*u21 + hn2*u22          ) +
00148                 anu*(-sint*(v11*t21 + v12*t22 + v13*t23) +
00149                     omcost*(v11*an1 + v12*an2 + v13*an3) -
00150               tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
00151   theJacobian(2,1) *= cosl1;
00152 
00153   theJacobian(2,2) = cost*(u11*u21 + u12*u22          ) +
00154                       sint*(hu1*u21 + hu2*u22          ) +
00155                     omcost*(hn1*u11 + hn2*u12          ) *
00156                            (hn1*u21 + hn2*u22          ) +
00157                 anu*(-sint*(u11*t21 + u12*t22          ) +
00158                     omcost*(u11*an1 + u12*an2          ) -
00159               tmsint*gamma*(hn1*u11 + hn2*u12          ) );
00160   theJacobian(2,2) *= cosl1*cosl0;
00161 
00162   theJacobian(2,3) = -q*anu*(u11*t21 + u12*t22          )*cosl1;
00163 
00164   theJacobian(2,4) = -q*anu*(v11*t21 + v12*t22 + v13*t23)*cosl1;
00165 
00166   //   yt
00167 
00168   //double cutCriterion = abs(s/fts.momentum().mag());
00169   double cutCriterion = fabs(s/globalParameters.momentum().mag());
00170   const double limit = 5.; // valid for propagations with effectively float precision
00171 
00172   if (cutCriterion > limit) {
00173     double pp = 1./qbp;
00174     theJacobian(3,0) = pp*(u21*dx1 + u22*dx2            );
00175     theJacobian(4,0) = pp*(v21*dx1 + v22*dx2 + v23*dx3);
00176   }
00177   else {
00178     double hp11 = hn2*t13 - hn3*t12;
00179     double hp12 = hn3*t11 - hn1*t13;
00180     double hp13 = hn1*t12 - hn2*t11;
00181     double temp1 = hp11*u21 + hp12*u22;
00182     double s2 = s*s;
00183     double secondOrder41 = 0.5 * qp * temp1 * s2;
00184     double ghnmp1 = gamma*hn1 - t11;
00185     double ghnmp2 = gamma*hn2 - t12;
00186     double ghnmp3 = gamma*hn3 - t13;
00187     double temp2 = ghnmp1*u21 + ghnmp2*u22;
00188     double s3 = s2 * s;
00189     double s4 = s3 * s;
00190     double h1 = h.mag();
00191     double h2 = h1 * h1;
00192     double h3 = h2 * h1;
00193     double qbp2 = qbp * qbp;
00194     double thirdOrder41 = 1./3 * h2 * s3 * qbp * temp2;
00195     double fourthOrder41 = 1./8 * h3 * s4 * qbp2 * temp1;
00196     theJacobian(3,0) = secondOrder41 + thirdOrder41 + fourthOrder41;
00197 
00198     double temp3 = hp11*v21 + hp12*v22 + hp13*v23;
00199     double secondOrder51 = 0.5 * qp * temp3 * s2;
00200     double temp4 = ghnmp1*v21 + ghnmp2*v22 + ghnmp3*v23;
00201     double thirdOrder51 = 1./3 * h2 * s3 * qbp * temp4;
00202     double fourthOrder51 = 1./8 * h3 * s4 * qbp2 * temp3;
00203     theJacobian(4,0) = secondOrder51 + thirdOrder51 + fourthOrder51;
00204   }
00205 
00206   theJacobian(3,1) = (sint*(v11*u21 + v12*u22          ) +
00207                      omcost*(hv1*u21 + hv2*u22          ) +
00208                      tmsint*(hn1*u21 + hn2*u22          ) *
00209                             (hn1*v11 + hn2*v12 + hn3*v13))/q;
00210 
00211   theJacobian(3,2) = (sint*(u11*u21 + u12*u22          ) +
00212                      omcost*(hu1*u21 + hu2*u22          ) +
00213                      tmsint*(hn1*u21 + hn2*u22          ) *
00214                             (hn1*u11 + hn2*u12          ))*cosl0/q;
00215 
00216   theJacobian(3,3) = (u11*u21 + u12*u22          );
00217   
00218   theJacobian(3,4) = (v11*u21 + v12*u22          );
00219 
00220   //   zt
00221 
00222   theJacobian(4,1) = (sint*(v11*v21 + v12*v22 + v13*v23) +
00223                      omcost*(hv1*v21 + hv2*v22 + hv3*v23) +
00224                      tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
00225                             (hn1*v11 + hn2*v12 + hn3*v13))/q;
00226 
00227   theJacobian(4,2) = (sint*(u11*v21 + u12*v22          ) +
00228                      omcost*(hu1*v21 + hu2*v22 + hu3*v23) +
00229                      tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
00230                             (hn1*u11 + hn2*u12          ))*cosl0/q;
00231 
00232   theJacobian(4,3) = (u11*v21 + u12*v22          );
00233 
00234   theJacobian(4,4) = (v11*v21 + v12*v22 + v13*v23);
00235   // end of TRPRFN
00236 }

void AnalyticalCurvilinearJacobian::computeStraightLineJacobian	(	const GlobalTrajectoryParameters &	globalParameters,
		const GlobalPoint &	x,
		const GlobalVector &	p,
		const double &	s
	)			`[private]`

straight line approximation

Definition at line 240 of file AnalyticalCurvilinearJacobian.cc.

References GlobalTrajectoryParameters::momentum(), p1, PV3DBase< T, PVType, FrameType >::perp(), and Vector3DBase< T, FrameTag >::unit().

00242 {
00243   //
00244   // matrix: elements =1 on diagonal and =0 are already set
00245   // in initialisation
00246   //
00247   GlobalVector p1 = globalParameters.momentum().unit();
00248   double cosl0 = p1.perp();
00249   theJacobian(3,2) = cosl0 * s;
00250   theJacobian(4,1) = s;
00251 }

virtual const AlgebraicMatrix55& AnalyticalCurvilinearJacobian::jacobian ( ) const [inline, virtual]

Implements CurvilinearJacobian.

Definition at line 38 of file AnalyticalCurvilinearJacobian.h.

References theJacobian.

00038 {return theJacobian;}

virtual const AlgebraicMatrix AnalyticalCurvilinearJacobian::jacobian_old ( ) const [inline, virtual]

Implements CurvilinearJacobian.

Definition at line 39 of file AnalyticalCurvilinearJacobian.h.

References asHepMatrix(), and theJacobian.

00039 {return asHepMatrix(theJacobian);}

AlgebraicMatrix55 AnalyticalCurvilinearJacobian::theJacobian [private]

Definition at line 52 of file AnalyticalCurvilinearJacobian.h.

Referenced by jacobian(), and jacobian_old().


Public Member Functions
	AnalyticalCurvilinearJacobian (const GlobalTrajectoryParameters &globalParameters, const GlobalPoint &x, const GlobalVector &p, const GlobalVector &theFieldInInverseGeV, const double &s)
	new: give Field as a parameter
	AnalyticalCurvilinearJacobian (const GlobalTrajectoryParameters &globalParameters, const GlobalPoint &x, const GlobalVector &p, const double &s)
	get Field at starting state (internally)
virtual const AlgebraicMatrix55 &	jacobian () const
virtual const AlgebraicMatrix	jacobian_old () const
virtual	~AnalyticalCurvilinearJacobian ()
Private Member Functions
void	computeFullJacobian (const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const GlobalVector &, const double &s)
	result for non-vanishing curvature
void	computeStraightLineJacobian (const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
	straight line approximation
Private Attributes
AlgebraicMatrix55	theJacobian

AnalyticalCurvilinearJacobian Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation