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:71
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:71
AnalyticalCurvilinearJacobian::computeStraightLineJacobian
void computeStraightLineJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
straight line approximation
Definition: AnalyticalCurvilinearJacobian.cc:311
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:58
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:71
AnalyticalCurvilinearJacobian::computeStraightLineJacobian
void computeStraightLineJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
straight line approximation
Definition: AnalyticalCurvilinearJacobian.cc:311
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:58

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 58 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(), x, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

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

254  {
255  /*
256  * origin TRPROP
257  *
258  C *** ERROR PROPAGATION ALONG A PARTICLE TRAJECTORY IN A MAGNETIC FIELD
259  C ROUTINE ASSUMES THAT IN THE INTERVAL (X1,X2) THE QUANTITIES 1/P
260  C AND (HX,HY,HZ) ARE RATHER CONSTANT. DELTA(PHI) MUST NOT BE TOO LARGE
261  C
262  C Authors: A. Haas and W. Wittek
263  C
264 
265  */
266 
267 
268  double qbp = globalParameters.signedInverseMomentum();
269  double absS = s;
270 
271  // average momentum
272  GlobalVector tn = (globalParameters.momentum()+p).unit();
273  double sinl = tn.z();
274  double cosl = std::sqrt(1.-sinl*sinl);
275  double cosl1 = 1./cosl;
276  double tgl=sinl*cosl1;
277  double sinp = tn.y()*cosl1;
278  double cosp = tn.x()*cosl1;
279 
280  // define average magnetic field and gradient
281  // at initial point - inlike TRPROP
282  double b0= h.x()*cosp+h.y()*sinp;
283  double b2=-h.x()*sinp+h.y()*cosp;
284  double b3=-b0*sinl+h.z()*cosl;
285 
286  theJacobian(3,2)=absS*cosl;
287  theJacobian(4,1)=absS;
288 
289 
290  theJacobian(1,0) = absS*b2;
291  //if ( qbp<0) theJacobian(1,0) = -theJacobian(1,0);
292  theJacobian(1,2) = -b0*(absS*qbp);
293  theJacobian(1,3) = b3*(b2*qbp*(absS*qbp));
294  theJacobian(1,4) = -b2*(b2*qbp*(absS*qbp));
295 
296  theJacobian(2,0) = -absS*b3*cosl1;
297  // if ( qbp<0) theJacobian(2,0) = -theJacobian(2,0);
298  theJacobian(2,1) = b0*(absS*qbp)*cosl1*cosl1;
299  theJacobian(2,2) = 1.+tgl*b2*(absS*qbp);
300  theJacobian(2,3) = -b3*(b3*qbp*(absS*qbp)*cosl1);
301  theJacobian(2,4) = b2*(b3*qbp*(absS*qbp)*cosl1);
302 
303  theJacobian(3,4) = -b3*tgl*(absS*qbp);
304  theJacobian(4,3) = b3*tgl*(absS*qbp);
305 
306 
307 }
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 311 of file AnalyticalCurvilinearJacobian.cc.

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

313 {
314  //
315  // matrix: elements =1 on diagonal and =0 are already set
316  // in initialisation
317  //
318  GlobalVector p1 = globalParameters.momentum().unit();
319  double cosl0 = p1.perp();
320  theJacobian(3,2) = cosl0 * s;
321  theJacobian(4,1) = s;
322 }
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
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 AnalyticalTrajectoryExtrapolatorToLine::extrapolateSingleState(), TSCBLBuilderNoMaterial::operator()(), 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