CMS 3D CMS Logo

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 AlgebraicMatrix55jacobian () 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() [1/3]

AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian ( )
inline

default constructor (for tests)

Definition at line 24 of file AnalyticalCurvilinearJacobian.h.

24 : theJacobian(ROOT::Math::SMatrixNoInit()) {}

◆ AnalyticalCurvilinearJacobian() [2/3]

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 computeFullJacobian(), computeStraightLineJacobian(), MillePedeFileConverter_cfg::e, h, GlobalTrajectoryParameters::magneticFieldInInverseGeV(), AlCaHLTBitMon_ParallelJobs::p, alignCSCRings::s, GlobalTrajectoryParameters::transverseCurvature(), and x.

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  GlobalVector h = globalParameters.magneticFieldInInverseGeV();
17  computeFullJacobian(globalParameters, x, p, h, s);
18  }
19  //
20  // straight line approximation, error in RPhi about 0.1um
21  //
22  else
23  computeStraightLineJacobian(globalParameters, x, p, s);
24  //dbg::dbg_trace(1,"ACJ1", globalParameters.vector(),x,p,s,theJacobian);
25 }
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
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4

◆ AnalyticalCurvilinearJacobian() [3/3]

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

References computeFullJacobian(), computeStraightLineJacobian(), MillePedeFileConverter_cfg::e, h, AlCaHLTBitMon_ParallelJobs::p, alignCSCRings::s, GlobalTrajectoryParameters::transverseCurvature(), and x.

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
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4

Member Function Documentation

◆ computeFullJacobian()

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 PVValHelper::dx, CustomPhysics_cfi::gamma, h, mps_fire::i, remoteMonitoring_LASER_era2018_cfg::limit, PV3DBase< T, PVType, FrameType >::mag(), GlobalTrajectoryParameters::momentum(), AlCaHLTBitMon_ParallelJobs::p, LaserDQM_cfg::p1, SiStripOfflineCRack_cfg::p2, GlobalTrajectoryParameters::position(), createTree::pp, submitPVResolutionJobs::q, alignCSCRings::s, photons_cff::s4, GlobalTrajectoryParameters::signedInverseMomentum(), mathSSE::sqrt(), theJacobian, theta(), Vector3DBase< T, FrameTag >::unit(), x, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by AnalyticalCurvilinearJacobian().

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

◆ computeInfinitesimalJacobian()

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 b0, b2, b3, h, GlobalTrajectoryParameters::momentum(), AlCaHLTBitMon_ParallelJobs::p, alignCSCRings::s, GlobalTrajectoryParameters::signedInverseMomentum(), mathSSE::sqrt(), theJacobian, 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  double qbp = globalParameters.signedInverseMomentum();
268  double absS = s;
269 
270  // average momentum
271  GlobalVector tn = (globalParameters.momentum() + p).unit();
272  double sinl = tn.z();
273  double cosl = std::sqrt(1. - sinl * sinl);
274  double cosl1 = 1. / cosl;
275  double tgl = sinl * cosl1;
276  double sinp = tn.y() * cosl1;
277  double cosp = tn.x() * cosl1;
278 
279  // define average magnetic field and gradient
280  // at initial point - inlike TRPROP
281  double b0 = h.x() * cosp + h.y() * sinp;
282  double b2 = -h.x() * sinp + h.y() * cosp;
283  double b3 = -b0 * sinl + h.z() * cosl;
284 
286 
287  theJacobian(3, 2) = absS * cosl;
288  theJacobian(4, 1) = absS;
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 }
weight_default_t b3[10]
Definition: b3.h:9
T z() const
Definition: PV3DBase.h:61
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
T sqrt(T t)
Definition: SSEVec.h:19
weight_default_t b2[10]
Definition: b2.h:9
Basic3DVector unit() const
static constexpr float b0
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4

◆ computeStraightLineJacobian()

void AnalyticalCurvilinearJacobian::computeStraightLineJacobian ( const GlobalTrajectoryParameters globalParameters,
const GlobalPoint ,
const GlobalVector ,
const double &  s 
)

straight line approximation

Definition at line 307 of file AnalyticalCurvilinearJacobian.cc.

References GlobalTrajectoryParameters::momentum(), LaserDQM_cfg::p1, alignCSCRings::s, theJacobian, and Vector3DBase< T, FrameTag >::unit().

Referenced by AnalyticalCurvilinearJacobian().

310  {
311  //
312  // matrix: elements =1 on diagonal and =0 are already set
313  // in initialisation
314  //
316  GlobalVector p1 = globalParameters.momentum().unit();
317  double cosl0 = p1.perp();
318  theJacobian(3, 2) = cosl0 * s;
319  theJacobian(4, 1) = s;
320 }
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Vector3DBase unit() const
Definition: Vector3DBase.h:54

◆ jacobian()

const AlgebraicMatrix55& AnalyticalCurvilinearJacobian::jacobian ( ) const
inline

Member Data Documentation

◆ theJacobian

AlgebraicMatrix55 AnalyticalCurvilinearJacobian::theJacobian
private