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>

Inheritance diagram for AnalyticalCurvilinearJacobian:
CurvilinearJacobian

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...
 
virtual const AlgebraicMatrix55jacobian () const
 
virtual const AlgebraicMatrix jacobian_old () const
 
virtual ~AnalyticalCurvilinearJacobian ()
 
- Public Member Functions inherited from CurvilinearJacobian
 CurvilinearJacobian ()
 
virtual ~CurvilinearJacobian ()
 

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.

ROOT::Math::SMatrixIdentity AlgebraicMatrixID
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 h, GlobalTrajectoryParameters::magneticFieldInInverseGeV(), GlobalTrajectoryParameters::position(), and GlobalTrajectoryParameters::transverseCurvature().

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 }
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
void computeStraightLineJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
straight line approximation
GlobalVector magneticFieldInInverseGeV(const GlobalPoint &x) const
void computeFullJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const GlobalVector &, const double &s)
result for non-vanishing curvature
string s
Definition: asciidump.py:422
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
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 GlobalTrajectoryParameters::transverseCurvature().

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 }
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
void computeStraightLineJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
straight line approximation
void computeFullJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const GlobalVector &, const double &s)
result for non-vanishing curvature
string s
Definition: asciidump.py:422
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
virtual AnalyticalCurvilinearJacobian::~AnalyticalCurvilinearJacobian ( )
inlinevirtual

Definition at line 38 of file AnalyticalCurvilinearJacobian.h.

38 {}

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, asciidump::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().

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  double thirdOrder41 = 1./3 * h2 * s3 * qbp * temp2;
199  double fourthOrder41 = 1./8 * h3 * s4 * qbp2 * temp1;
200  theJacobian(3,0) = secondOrder41 + (thirdOrder41 + fourthOrder41);
201 
202  double temp3 = hp11*v21 + hp12*v22 + hp13*v23;
203  double secondOrder51 = 0.5 * qp * temp3 * s2;
204  double temp4 = ghnmp1*v21 + ghnmp2*v22 + ghnmp3*v23;
205  double thirdOrder51 = 1./3 * h2 * s3 * qbp * temp4;
206  double fourthOrder51 = 1./8 * h3 * s4 * qbp2 * temp3;
207  theJacobian(4,0) = secondOrder51 + (thirdOrder51 + fourthOrder51);
208  }
209 
210  theJacobian(3,1) = (sint*(v11*u21 + v12*u22 ) +
211  omcost*(hv1*u21 + hv2*u22 ) +
212  tmsint*(hn1*u21 + hn2*u22 ) *
213  (hn1*v11 + hn2*v12 + hn3*v13))/q;
214 
215  theJacobian(3,2) = (sint*(u11*u21 + u12*u22 ) +
216  omcost*(hu1*u21 + hu2*u22 ) +
217  tmsint*(hn1*u21 + hn2*u22 ) *
218  (hn1*u11 + hn2*u12 ))*cosl0/q;
219 
220  theJacobian(3,3) = (u11*u21 + u12*u22 );
221 
222  theJacobian(3,4) = (v11*u21 + v12*u22 );
223 
224  // zt
225 
226  theJacobian(4,1) = (sint*(v11*v21 + v12*v22 + v13*v23) +
227  omcost*(hv1*v21 + hv2*v22 + hv3*v23) +
228  tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
229  (hn1*v11 + hn2*v12 + hn3*v13))/q;
230 
231  theJacobian(4,2) = (sint*(u11*v21 + u12*v22 ) +
232  omcost*(hu1*v21 + hu2*v22 + hu3*v23) +
233  tmsint*(hn1*v21 + hn2*v22 + hn3*v23) *
234  (hn1*u11 + hn2*u12 ))*cosl0/q;
235 
236  theJacobian(4,3) = (u11*v21 + u12*v22 );
237 
238  theJacobian(4,4) = (v11*v21 + v12*v22 + v13*v23);
239  // end of TRPRFN
240 }
T perp() const
Definition: PV3DBase.h:66
tuple pp
Definition: createTree.py:15
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Theta< T > theta() const
T y() const
Definition: PV3DBase.h:57
tuple s2
Definition: indexGen.py:106
T mag() const
Definition: PV3DBase.h:61
T sqrt(T t)
Definition: SSEVec.h:28
T z() const
Definition: PV3DBase.h:58
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
double p2[4]
Definition: TauolaWrapper.h:90
Vector3DBase unit() const
Definition: Vector3DBase.h:57
double p1[4]
Definition: TauolaWrapper.h:89
string s
Definition: asciidump.py:422
Definition: DDAxes.h:10
T x() const
Definition: PV3DBase.h:56
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 245 of file AnalyticalCurvilinearJacobian.cc.

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

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

straight line approximation

Definition at line 306 of file AnalyticalCurvilinearJacobian.cc.

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

308 {
309  //
310  // matrix: elements =1 on diagonal and =0 are already set
311  // in initialisation
312  //
313  GlobalVector p1 = globalParameters.momentum().unit();
314  double cosl0 = p1.perp();
315  theJacobian(3,2) = cosl0 * s;
316  theJacobian(4,1) = s;
317 }
T perp() const
Definition: PV3DBase.h:66
Vector3DBase unit() const
Definition: Vector3DBase.h:57
double p1[4]
Definition: TauolaWrapper.h:89
string s
Definition: asciidump.py:422
virtual const AlgebraicMatrix55& AnalyticalCurvilinearJacobian::jacobian ( ) const
inlinevirtual
virtual const AlgebraicMatrix AnalyticalCurvilinearJacobian::jacobian_old ( ) const
inlinevirtual

Implements CurvilinearJacobian.

Definition at line 41 of file AnalyticalCurvilinearJacobian.h.

References asHepMatrix(), and theJacobian.

41 {return asHepMatrix(theJacobian);}
CLHEP::HepMatrix asHepMatrix(const ROOT::Math::SMatrix< double, N1, N2, typename ROOT::Math::MatRepStd< double, N1, N2 > > &rm)
Definition: Migration.h:49

Member Data Documentation

AlgebraicMatrix55 AnalyticalCurvilinearJacobian::theJacobian
private

Definition at line 58 of file AnalyticalCurvilinearJacobian.h.

Referenced by jacobian(), and jacobian_old().