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(ROOT::Math::SMatrixNoInit()) {}
AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian ( const GlobalTrajectoryParameters &  globalParameters,
const GlobalPoint &  x,
const GlobalVector &  p,
const double &  s 
)

get Field at starting state (internally)

Definition at line 7 of file AnalyticalCurvilinearJacobian.cc.

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

10  : theJacobian(AlgebraicMatrixID())
11 {
12  //
13  // helix: calculate full jacobian
14  //
15  if ( s*s*fabs(globalParameters.transverseCurvature())>1.e-5 ) {
16  // GlobalPoint xStart = globalParameters.position();
17  // GlobalVector h = globalParameters.magneticFieldInInverseGeV(xStart);
18  GlobalVector h = globalParameters.magneticFieldInInverseGeV();
19  computeFullJacobian(globalParameters,x,p,h,s);
20  }
21  //
22  // straight line approximation, error in RPhi about 0.1um
23  //
24  else
25  computeStraightLineJacobian(globalParameters,x,p,s);
26  //dbg::dbg_trace(1,"ACJ1", globalParameters.vector(),x,p,s,theJacobian);
27 }
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:80
AnalyticalCurvilinearJacobian::computeStraightLineJacobian
void computeStraightLineJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
straight line approximation
Definition: AnalyticalCurvilinearJacobian.cc:316
GlobalTrajectoryParameters::magneticFieldInInverseGeV
GlobalVector magneticFieldInInverseGeV(const GlobalPoint &x) const
Definition: GlobalTrajectoryParameters.cc:64
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:61
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 31 of file AnalyticalCurvilinearJacobian.cc.

References alignCSCRings::e, and GlobalTrajectoryParameters::transverseCurvature().

35  : theJacobian(AlgebraicMatrixID())
36 {
37  //
38  // helix: calculate full jacobian
39  //
40  if ( s*s*fabs(globalParameters.transverseCurvature())>1.e-5 )
41  computeFullJacobian(globalParameters,x,p,h,s);
42  //
43  // straight line approximation, error in RPhi about 0.1um
44  //
45  else
46  computeStraightLineJacobian(globalParameters,x,p,s);
47 
48  //dbg::dbg_trace(1,"ACJ2", globalParameters.vector(),x,p,s,theJacobian);
49 }
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:80
AnalyticalCurvilinearJacobian::computeStraightLineJacobian
void computeStraightLineJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
straight line approximation
Definition: AnalyticalCurvilinearJacobian.cc:316
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:61

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 61 of file AnalyticalCurvilinearJacobian.cc.

References prof2calltree::cost, i, 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(), mathSSE::sqrt(), theta(), Vector3DBase< T, FrameTag >::unit(), x, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

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

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

straight line approximation

Definition at line 316 of file AnalyticalCurvilinearJacobian.cc.

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

318 {
319  //
320  // matrix: elements =1 on diagonal and =0 are already set
321  // in initialisation
322  //
323  theJacobian = AlgebraicMatrixID();
324  GlobalVector p1 = globalParameters.momentum().unit();
325  double cosl0 = p1.perp();
326  theJacobian(3,2) = cosl0 * s;
327  theJacobian(4,1) = s;
328 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:80
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 analyticalErrorPropagation(), AnalyticalTrajectoryExtrapolatorToLine::extrapolateSingleState(), TSCBLBuilderNoMaterial::operator()(), TwoBodyDecayTrajectoryState::propagateSingleState(), and TrackKinematicStatePropagator::propagateToTheTransversePCACharged().

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