CMS 3D CMS Logo
AnalyticalCurvilinearJacobian.cc
Go to the documentation of this file.
1 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
2 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
3 #include <vdt/vdtMath.h>
4 
5 AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian(const GlobalTrajectoryParameters& globalParameters,
6  const GlobalPoint& x,
7  const GlobalVector& p,
8  const double& s)
9  : theJacobian(AlgebraicMatrixID()) {
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 }
26 
27 AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian(
28  const GlobalTrajectoryParameters& globalParameters,
29  const GlobalPoint& x,
30  const GlobalVector& p,
31  const GlobalVector& h, // h is the magnetic Field in Inverse GeV
32  const double& s)
33  : theJacobian(AlgebraicMatrixID()) {
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 }
47 
48 #if defined(USE_SSEVECT) && !defined(TRPRFN_SCALAR)
49 #include "AnalyticalCurvilinearJacobianSSE.icc"
50 #elif defined(USE_EXTVECT) && !defined(TRPRFN_SCALAR)
51 #include "AnalyticalCurvilinearJacobianEXT.icc"
52 
53 #else
54 
55 void AnalyticalCurvilinearJacobian::computeFullJacobian(const GlobalTrajectoryParameters& globalParameters,
56  const GlobalPoint& x,
57  const GlobalVector& p,
58  const GlobalVector& h,
59  const double& s) {
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 }
247 
248 #endif
249 
250 void AnalyticalCurvilinearJacobian::computeInfinitesimalJacobian(const GlobalTrajectoryParameters& globalParameters,
251  const GlobalPoint&,
252  const GlobalVector& p,
253  const GlobalVector& h,
254  const double& s) {
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 
285  theJacobian = AlgebraicMatrixID();
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 }
306 
307 void AnalyticalCurvilinearJacobian::computeStraightLineJacobian(const GlobalTrajectoryParameters& globalParameters,
308  const GlobalPoint&,
309  const GlobalVector&,
310  const double& s) {
311  //
312  // matrix: elements =1 on diagonal and =0 are already set
313  // in initialisation
314  //
315  theJacobian = AlgebraicMatrixID();
316  GlobalVector p1 = globalParameters.momentum().unit();
317  double cosl0 = p1.perp();
318  theJacobian(3, 2) = cosl0 * s;
319  theJacobian(4, 1) = s;
320 }
Vector3DBase
Definition: Vector3DBase.h:8
AnalyticalCurvilinearJacobian::AnalyticalCurvilinearJacobian
AnalyticalCurvilinearJacobian()
default constructor (for tests)
Definition: AnalyticalCurvilinearJacobian.h:24
mps_fire.i
i
Definition: mps_fire.py:428
GlobalTrajectoryParameters.h
GlobalTrajectoryParameters::position
GlobalPoint position() const
Definition: GlobalTrajectoryParameters.h:60
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
AnalyticalCurvilinearJacobian::computeInfinitesimalJacobian
void computeInfinitesimalJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const GlobalVector &, const double &s)
result for non-vanishing curvature and "small" step
Definition: AnalyticalCurvilinearJacobian.cc:250
CustomPhysics_cfi.gamma
gamma
Definition: CustomPhysics_cfi.py:17
AnalyticalCurvilinearJacobian.h
DDAxes::x
indexGen.s2
s2
Definition: indexGen.py:107
b2
static constexpr float b2
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:72
alignCSCRings.s
s
Definition: alignCSCRings.py:92
Vector3DBase::unit
Vector3DBase unit() const
Definition: Vector3DBase.h:54
AnalyticalCurvilinearJacobian::theJacobian
AlgebraicMatrix55 theJacobian
Definition: AnalyticalCurvilinearJacobian.h:54
AnalyticalCurvilinearJacobian::computeStraightLineJacobian
void computeStraightLineJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const double &s)
straight line approximation
Definition: AnalyticalCurvilinearJacobian.cc:307
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
p2
double p2[4]
Definition: TauolaWrapper.h:90
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
GlobalTrajectoryParameters
Definition: GlobalTrajectoryParameters.h:15
h
Point3DBase< float, GlobalTag >
GlobalTrajectoryParameters::momentum
GlobalVector momentum() const
Definition: GlobalTrajectoryParameters.h:65
GlobalTrajectoryParameters::transverseCurvature
float transverseCurvature() const
Definition: GlobalTrajectoryParameters.h:89
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
p1
double p1[4]
Definition: TauolaWrapper.h:89
submitPVResolutionJobs.q
q
Definition: submitPVResolutionJobs.py:84
b0
static constexpr float b0
Definition: L1EGammaCrystalsEmulatorProducer.cc:83
unit
Basic3DVector unit() const
Definition: Basic3DVectorLD.h:162
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:64
remoteMonitoring_LED_IterMethod_cfg.limit
limit
Definition: remoteMonitoring_LED_IterMethod_cfg.py:427
testProducerWithPsetDescEmpty_cfi.b3
b3
Definition: testProducerWithPsetDescEmpty_cfi.py:36
genVertex_cff.x
x
Definition: genVertex_cff.py:12
GlobalTrajectoryParameters::magneticFieldInInverseGeV
GlobalVector magneticFieldInInverseGeV(const GlobalPoint &x) const
Definition: GlobalTrajectoryParameters.cc:35
AnalyticalCurvilinearJacobian::computeFullJacobian
void computeFullJacobian(const GlobalTrajectoryParameters &, const GlobalPoint &, const GlobalVector &, const GlobalVector &, const double &s)
result for non-vanishing curvature
Definition: AnalyticalCurvilinearJacobian.cc:55
createTree.pp
pp
Definition: createTree.py:17
GlobalTrajectoryParameters::signedInverseMomentum
float signedInverseMomentum() const
Definition: GlobalTrajectoryParameters.h:77
PVValHelper::dx
Definition: PVValidationHelpers.h:49
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16