CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Attributes
AnalyticalCurvilinearJacobian Class Reference

#include <AnalyticalCurvilinearJacobian.h>

Public Member Functions

 AnalyticalCurvilinearJacobian ()
 default constructor (for tests) More...
 
 AnalyticalCurvilinearJacobian (const GlobalTrajectoryParameters &globalParameters, const GlobalPoint &x, const GlobalVector &p, const double &s)
 get Field at starting state (internally) More...
 
 AnalyticalCurvilinearJacobian (const GlobalTrajectoryParameters &globalParameters, const GlobalPoint &x, const GlobalVector &p, const GlobalVector &theFieldInInverseGeV, const double &s)
 new: give Field as a parameter More...
 
void computeFullJacobian (const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const GlobalVector &, const double &s)
 result for non-vanishing curvature More...
 
void computeInfinitesimalJacobian (const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const GlobalVector &, const double &s)
 result for non-vanishing curvature and "small" step More...
 
void computeStraightLineJacobian (const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
 straight line approximation More...
 
const AlgebraicMatrix55 & jacobian () const
 

Private Attributes

AlgebraicMatrix55 theJacobian
 

Detailed Description

Creating Jacobian of transformation within the curvilinear frame. The basic functionality of this class is to provide the (analytical) Jacobian matrix of the transformation within the curvilinear frame from the state defined by globalParameters to the state defined by x and p. This Jacobian can then be used to yield the corresponding error propagation. The current implementation is based on the original derivations by W. Wittek. However, due to the implicit float precision, two terms ((4,1) and (5,1)) have been modified in order to make the calculations more stable in a numerical sense.

Definition at line 21 of file AnalyticalCurvilinearJacobian.h.

Constructor & Destructor Documentation

AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian ( )
inline

default constructor (for tests)

Definition at line 24 of file AnalyticalCurvilinearJacobian.h.

24 : theJacobian(AlgebraicMatrixID()){}
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:70
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 alignCSCRings::e, h, GlobalTrajectoryParameters::magneticFieldInInverseGeV(), GlobalTrajectoryParameters::position(), and GlobalTrajectoryParameters::transverseCurvature().

8  : theJacobian(AlgebraicMatrixID())
9 {
10  //
11  // helix: calculate full jacobian
12  //
13  if ( s*s*fabs(globalParameters.transverseCurvature())>1.e-5 ) {
14  GlobalPoint xStart = globalParameters.position();
15  GlobalVector h = globalParameters.magneticFieldInInverseGeV(xStart);
16  computeFullJacobian(globalParameters,x,p,h,s);
17  }
18  //
19  // straight line approximation, error in RPhi about 0.1um
20  //
21  else
22  computeStraightLineJacobian(globalParameters,x,p,s);
23  //dbg::dbg_trace(1,"ACJ1", globalParameters.vector(),x,p,s,theJacobian);
24 }
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:70
AnalyticalCurvilinearJacobian::computeStraightLineJacobian
void computeStraightLineJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
straight line approximation
Definition: AnalyticalCurvilinearJacobian.cc:308
GlobalTrajectoryParameters::magneticFieldInInverseGeV
GlobalVector magneticFieldInInverseGeV(const GlobalPoint &x) const
Definition: GlobalTrajectoryParameters.cc:28
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
AnalyticalCurvilinearJacobian::computeFullJacobian
void computeFullJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const GlobalVector &, const double &s)
result for non-vanishing curvature
Definition: AnalyticalCurvilinearJacobian.cc:55
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 alignCSCRings::e, and GlobalTrajectoryParameters::transverseCurvature().

32  : theJacobian(AlgebraicMatrixID())
33 {
34  //
35  // helix: calculate full jacobian
36  //
37  if ( s*s*fabs(globalParameters.transverseCurvature())>1.e-5 )
38  computeFullJacobian(globalParameters,x,p,h,s);
39  //
40  // straight line approximation, error in RPhi about 0.1um
41  //
42  else
43  computeStraightLineJacobian(globalParameters,x,p,s);
44 
45  //dbg::dbg_trace(1,"ACJ2", globalParameters.vector(),x,p,s,theJacobian);
46 }
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:70
AnalyticalCurvilinearJacobian::computeStraightLineJacobian
void computeStraightLineJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
straight line approximation
Definition: AnalyticalCurvilinearJacobian.cc:308
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
AnalyticalCurvilinearJacobian::computeFullJacobian
void computeFullJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const GlobalVector &, const double &s)
result for non-vanishing curvature
Definition: AnalyticalCurvilinearJacobian.cc:55

Member Function Documentation

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

result for non-vanishing curvature

Definition at line 55 of file AnalyticalCurvilinearJacobian.cc.

References funct::cos(), prof2calltree::cost, MessageLogger_cff::limit, PV3DBase< T, PVType, FrameType >::mag(), GlobalTrajectoryParameters::momentum(), p1, p2, PV3DBase< T, PVType, FrameType >::perp(), GlobalTrajectoryParameters::position(), createTree::pp, lumiQueryAPI::q, alignCSCRings::s, indexGen::s2, GlobalTrajectoryParameters::signedInverseMomentum(), funct::sin(), mathSSE::sqrt(), theta(), Vector3DBase< T, FrameTag >::unit(), PV3DBase< T, PVType, FrameType >::x(), vdt::x, PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

60 {
61  //GlobalVector p1 = fts.momentum().unit();
62  GlobalVector p1 = globalParameters.momentum().unit();
63  GlobalVector p2 = p.unit();
64  //GlobalPoint xStart = fts.position();
65  GlobalPoint xStart = globalParameters.position();
66  GlobalVector dx = xStart - x;
67  //GlobalVector h = MagneticField::inInverseGeV(xStart);
68  // Martijn: field is now given as parameter.. GlobalVector h = globalParameters.magneticFieldInInverseGeV(xStart);
69 
70  //double qbp = fts.signedInverseMomentum();
71  double qbp = globalParameters.signedInverseMomentum();
72  double absS = s;
73 
74  // calculate transport matrix
75  // Origin: TRPRFN
76  double t11 = p1.x(); double t12 = p1.y(); double t13 = p1.z();
77  double t21 = p2.x(); double t22 = p2.y(); double t23 = p2.z();
78  double cosl0 = p1.perp(); double cosl1 = 1./p2.perp();
79  //AlgebraicMatrix a(5,5,1);
80  // define average magnetic field and gradient
81  // at initial point - inlike TRPRFN
82  GlobalVector hn = h.unit();
83  double qp = -h.mag();
84 // double q = -h.mag()*qbp;
85  double q = qp*qbp;
86  double theta = q*absS; double sint = sin(theta); double cost = cos(theta);
87  double hn1 = hn.x(); double hn2 = hn.y(); double hn3 = hn.z();
88  double dx1 = dx.x(); double dx2 = dx.y(); double dx3 = dx.z();
89  double gamma = hn1*t21 + hn2*t22 + hn3*t23;
90  double an1 = hn2*t23 - hn3*t22;
91  double an2 = hn3*t21 - hn1*t23;
92  double an3 = hn1*t22 - hn2*t21;
93  double au = 1./sqrt(t11*t11 + t12*t12);
94  double u11 = -au*t12; double u12 = au*t11;
95  double v11 = -t13*u12; double v12 = t13*u11; double v13 = t11*u12 - t12*u11;
96  au = 1./sqrt(t21*t21 + t22*t22);
97  double u21 = -au*t22; double u22 = au*t21;
98  double v21 = -t23*u22; double v22 = t23*u21; double v23 = t21*u22 - t22*u21;
99  // now prepare the transport matrix
100  // pp only needed in high-p case (WA)
101 // double pp = 1./qbp;
103 // moved up (where -h.mag() is needed()
104 // double qp = q*pp;
105  double anv = -(hn1*u21 + hn2*u22 );
106  double anu = (hn1*v21 + hn2*v22 + hn3*v23);
107  double omcost = 1. - cost; double tmsint = theta - sint;
108 
109  double hu1 = - hn3*u12;
110  double hu2 = hn3*u11;
111  double hu3 = hn1*u12 - hn2*u11;
112 
113  double hv1 = hn2*v13 - hn3*v12;
114  double hv2 = hn3*v11 - hn1*v13;
115  double hv3 = hn1*v12 - hn2*v11;
116 
117  // 1/p - doesn't change since |p1| = |p2|
118 
119  // lambda
120 
121  theJacobian(1,0) = -qp*anv*(t21*dx1 + t22*dx2 + t23*dx3);
122 
123  theJacobian(1,1) = cost*(v11*v21 + v12*v22 + v13*v23) +
124  sint*(hv1*v21 + hv2*v22 + hv3*v23) +
125  omcost*(hn1*v11 + hn2*v12 + hn3*v13) *
126  (hn1*v21 + hn2*v22 + hn3*v23) +
127  anv*(-sint*(v11*t21 + v12*t22 + v13*t23) +
128  omcost*(v11*an1 + v12*an2 + v13*an3) -
129  tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
130 
131  theJacobian(1,2) = cost*(u11*v21 + u12*v22 ) +
132  sint*(hu1*v21 + hu2*v22 + hu3*v23) +
133  omcost*(hn1*u11 + hn2*u12 ) *
134  (hn1*v21 + hn2*v22 + hn3*v23) +
135  anv*(-sint*(u11*t21 + u12*t22 ) +
136  omcost*(u11*an1 + u12*an2 ) -
137  tmsint*gamma*(hn1*u11 + hn2*u12 ) );
138  theJacobian(1,2) *= cosl0;
139 
140  theJacobian(1,3) = -q*anv*(u11*t21 + u12*t22 );
141 
142  theJacobian(1,4) = -q*anv*(v11*t21 + v12*t22 + v13*t23);
143 
144  // phi
145 
146  theJacobian(2,0) = -qp*anu*(t21*dx1 + t22*dx2 + t23*dx3)*cosl1;
147 
148  theJacobian(2,1) = cost*(v11*u21 + v12*u22 ) +
149  sint*(hv1*u21 + hv2*u22 ) +
150  omcost*(hn1*v11 + hn2*v12 + hn3*v13) *
151  (hn1*u21 + hn2*u22 ) +
152  anu*(-sint*(v11*t21 + v12*t22 + v13*t23) +
153  omcost*(v11*an1 + v12*an2 + v13*an3) -
154  tmsint*gamma*(hn1*v11 + hn2*v12 + hn3*v13) );
155  theJacobian(2,1) *= cosl1;
156 
157  theJacobian(2,2) = cost*(u11*u21 + u12*u22 ) +
158  sint*(hu1*u21 + hu2*u22 ) +
159  omcost*(hn1*u11 + hn2*u12 ) *
160  (hn1*u21 + hn2*u22 ) +
161  anu*(-sint*(u11*t21 + u12*t22 ) +
162  omcost*(u11*an1 + u12*an2 ) -
163  tmsint*gamma*(hn1*u11 + hn2*u12 ) );
164  theJacobian(2,2) *= cosl1*cosl0;
165 
166  theJacobian(2,3) = -q*anu*(u11*t21 + u12*t22 )*cosl1;
167 
168  theJacobian(2,4) = -q*anu*(v11*t21 + v12*t22 + v13*t23)*cosl1;
169 
170  // yt
171 
172  //double cutCriterion = abs(s/fts.momentum().mag());
173  double cutCriterion = fabs(s/globalParameters.momentum().mag());
174  const double limit = 5.; // valid for propagations with effectively float precision
175 
176  if (cutCriterion > limit) {
177  double pp = 1./qbp;
178  theJacobian(3,0) = pp*(u21*dx1 + u22*dx2 );
179  theJacobian(4,0) = pp*(v21*dx1 + v22*dx2 + v23*dx3);
180  }
181  else {
182  double hp11 = hn2*t13 - hn3*t12;
183  double hp12 = hn3*t11 - hn1*t13;
184  double hp13 = hn1*t12 - hn2*t11;
185  double temp1 = hp11*u21 + hp12*u22;
186  double s2 = s*s;
187  double secondOrder41 = 0.5 * qp * temp1 * s2;
188  double ghnmp1 = gamma*hn1 - t11;
189  double ghnmp2 = gamma*hn2 - t12;
190  double ghnmp3 = gamma*hn3 - t13;
191  double temp2 = ghnmp1*u21 + ghnmp2*u22;
192  double s3 = s2 * s;
193  double s4 = s3 * s;
194  double h1 = h.mag();
195  double h2 = h1 * h1;
196  double h3 = h2 * h1;
197  double qbp2 = qbp * qbp;
198  // s*qp*s* (qp*s *qbp)
199  double thirdOrder41 = 1./3 * h2 * s3 * qbp * temp2;
200  // -qp * s * qbp * above
201  double fourthOrder41 = 1./8 * h3 * s4 * qbp2 * temp1;
202  theJacobian(3,0) = secondOrder41 + (thirdOrder41 + fourthOrder41);
203 
204  double temp3 = hp11*v21 + hp12*v22 + hp13*v23;
205  double secondOrder51 = 0.5 * qp * temp3 * s2;
206  double temp4 = ghnmp1*v21 + ghnmp2*v22 + ghnmp3*v23;
207  double thirdOrder51 = 1./3 * h2 * s3 * qbp * temp4;
208  double fourthOrder51 = 1./8 * h3 * s4 * qbp2 * temp3;
209  theJacobian(4,0) = secondOrder51 + (thirdOrder51 + fourthOrder51);
210  }
211 
212  theJacobian(3,1) = (sint*(v11*u21 + v12*u22 ) +
213  omcost*(hv1*u21 + hv2*u22 ) +
214  tmsint*(hn1*u21 + hn2*u22 ) *
215  (hn1*v11 + hn2*v12 + hn3*v13))/q;
216 
217  theJacobian(3,2) = (sint*(u11*u21 + u12*u22 ) +
218  omcost*(hu1*u21 + hu2*u22 ) +
219  tmsint*(hn1*u21 + hn2*u22 ) *
220  (hn1*u11 + hn2*u12 ))*cosl0/q;
221 
222  theJacobian(3,3) = (u11*u21 + u12*u22 );
223 
224  theJacobian(3,4) = (v11*u21 + v12*u22 );
225 
226  // zt
227 
228  theJacobian(4,1) = (sint*(v11*v21 + v12*v22 + v13*v23) +
229  omcost*(hv1*v21 + hv2*v22 + hv3*v23) +
230  tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
231  (hn1*v11 + hn2*v12 + hn3*v13))/q;
232 
233  theJacobian(4,2) = (sint*(u11*v21 + u12*v22 ) +
234  omcost*(hu1*v21 + hu2*v22 + hu3*v23) +
235  tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
236  (hn1*u11 + hn2*u12 ))*cosl0/q;
237 
238  theJacobian(4,3) = (u11*v21 + u12*v22 );
239 
240  theJacobian(4,4) = (v11*v21 + v12*v22 + v13*v23);
241  // end of TRPRFN
242 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:71
createTree.pp
tuple pp
Definition: createTree.py:15
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:165
PV3DBase::y
T y() const
Definition: PV3DBase.h:62
MessageLogger_cff.limit
tuple limit
Definition: MessageLogger_cff.py:11
indexGen.s2
tuple s2
Definition: indexGen.py:106
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:66
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
PV3DBase::z
T z() const
Definition: PV3DBase.h:63
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
p2
double p2[4]
Definition: TauolaWrapper.h:90
Vector3DBase::unit
Vector3DBase unit() const
Definition: Vector3DBase.h:57
p1
double p1[4]
Definition: TauolaWrapper.h:89
vdt::x
x
Definition: VDTMath.h:216
PV3DBase::x
T x() const
Definition: PV3DBase.h:61
void AnalyticalCurvilinearJacobian::computeInfinitesimalJacobian ( const GlobalTrajectoryParameters &  globalParameters,
const GlobalPoint &  ,
const GlobalVector &  p,
const GlobalVector &  h,
const double &  s 
)

result for non-vanishing curvature and "small" step

Definition at line 247 of file AnalyticalCurvilinearJacobian.cc.

References GlobalTrajectoryParameters::momentum(), AlCaHLTBitMon_ParallelJobs::p, alignCSCRings::s, GlobalTrajectoryParameters::signedInverseMomentum(), mathSSE::sqrt(), csvLumiCalc::unit, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

251  {
252  /*
253  * origin TRPROP
254  *
255  C *** ERROR PROPAGATION ALONG A PARTICLE TRAJECTORY IN A MAGNETIC FIELD
256  C ROUTINE ASSUMES THAT IN THE INTERVAL (X1,X2) THE QUANTITIES 1/P
257  C AND (HX,HY,HZ) ARE RATHER CONSTANT. DELTA(PHI) MUST NOT BE TOO LARGE
258  C
259  C Authors: A. Haas and W. Wittek
260  C
261 
262  */
263 
264 
265  double qbp = globalParameters.signedInverseMomentum();
266  double absS = s;
267 
268  // average momentum
269  GlobalVector tn = (globalParameters.momentum()+p).unit();
270  double sinl = tn.z();
271  double cosl = std::sqrt(1.-sinl*sinl);
272  double cosl1 = 1./cosl;
273  double tgl=sinl*cosl1;
274  double sinp = tn.y()*cosl1;
275  double cosp = tn.x()*cosl1;
276 
277  // define average magnetic field and gradient
278  // at initial point - inlike TRPROP
279  double b0= h.x()*cosp+h.y()*sinp;
280  double b2=-h.x()*sinp+h.y()*cosp;
281  double b3=-b0*sinl+h.z()*cosl;
282 
283  theJacobian(3,2)=absS*cosl;
284  theJacobian(4,1)=absS;
285 
286 
287  theJacobian(1,0) = absS*b2;
288  //if ( qbp<0) theJacobian(1,0) = -theJacobian(1,0);
289  theJacobian(1,2) = -b0*(absS*qbp);
290  theJacobian(1,3) = b3*(b2*qbp*(absS*qbp));
291  theJacobian(1,4) = -b2*(b2*qbp*(absS*qbp));
292 
293  theJacobian(2,0) = -absS*b3*cosl1;
294  // if ( qbp<0) theJacobian(2,0) = -theJacobian(2,0);
295  theJacobian(2,1) = b0*(absS*qbp)*cosl1*cosl1;
296  theJacobian(2,2) = 1.+tgl*b2*(absS*qbp);
297  theJacobian(2,3) = -b3*(b3*qbp*(absS*qbp)*cosl1);
298  theJacobian(2,4) = b2*(b3*qbp*(absS*qbp)*cosl1);
299 
300  theJacobian(3,4) = -b3*tgl*(absS*qbp);
301  theJacobian(4,3) = b3*tgl*(absS*qbp);
302 
303 
304 }
PV3DBase::y
T y() const
Definition: PV3DBase.h:62
csvLumiCalc.unit
string unit
Definition: csvLumiCalc.py:46
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
PV3DBase::z
T z() const
Definition: PV3DBase.h:63
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
PV3DBase::x
T x() const
Definition: PV3DBase.h:61
void AnalyticalCurvilinearJacobian::computeStraightLineJacobian ( const GlobalTrajectoryParameters &  globalParameters,
const GlobalPoint &  x,
const GlobalVector &  p,
const double &  s 
)

straight line approximation

Definition at line 308 of file AnalyticalCurvilinearJacobian.cc.

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

310 {
311  //
312  // matrix: elements =1 on diagonal and =0 are already set
313  // in initialisation
314  //
315  GlobalVector p1 = globalParameters.momentum().unit();
316  double cosl0 = p1.perp();
317  theJacobian(3,2) = cosl0 * s;
318  theJacobian(4,1) = s;
319 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:71
Vector3DBase::unit
Vector3DBase unit() const
Definition: Vector3DBase.h:57
p1
double p1[4]
Definition: TauolaWrapper.h:89
const AlgebraicMatrix55& AnalyticalCurvilinearJacobian::jacobian ( ) const
inline

Definition at line 55 of file AnalyticalCurvilinearJacobian.h.

References theJacobian.

Referenced by TSCPBuilderNoMaterial::createFTSatTransverseImpactPointCharged(), TSCPBuilderNoMaterial::createFTSatTransverseImpactPointNeutral(), AnalyticalTrajectoryExtrapolatorToLine::extrapolateSingleState(), TSCBLBuilderNoMaterial::operator()(), AnalyticalPropagator::propagatedStateWithPath(), and TwoBodyDecayTrajectoryState::propagateSingleState().

55 {return theJacobian;}

Member Data Documentation

AlgebraicMatrix55 AnalyticalCurvilinearJacobian::theJacobian
private

Definition at line 59 of file AnalyticalCurvilinearJacobian.h.

Referenced by jacobian().

Generated for CMSSW Reference Manual by doxygen 1.8.5