CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
GhostTrackPrediction.cc
Go to the documentation of this file.
1 #include <cmath>
2 
3 #include <Math/SMatrix.h>
4 #include <Math/MatrixFunctions.h>
5 
10 
15 
17 
18 using namespace reco;
19 
20 namespace {
21  static inline double sqr(double arg) { return arg * arg; }
22 
23  using namespace ROOT::Math;
24 
25  typedef SMatrix<double, 3, 4> Matrix34;
26  typedef SMatrix<double, 4, 5> Matrix45;
27  typedef SMatrix<double, 5, 4> Matrix54;
28  typedef SMatrix<double, 4, 6> Matrix46;
29  typedef SMatrix<double, 6, 4> Matrix64;
30  typedef SMatrix<double, 6, 6, MatRepSym<double, 6> > Matrix6S;
31 }
32 
34  const CurvilinearTrajectoryParameters &trajectory)
35 {
37  trajectory.yT() / std::cos(trajectory.lambda()),
38  trajectory.xT(),
39  std::tan(trajectory.lambda()),
40  trajectory.phi());
41 }
42 
44  const GlobalTrajectoryParameters &trajectory)
45 {
47  trajectory.position(), trajectory.momentum(), 0));
48 }
49 
51  const GhostTrackPrediction::Vector &pred,
53 {
54  using namespace ROOT::Math;
55 
56  double rho2 = pred[2] * pred[2] + 1.;
57  double rho = std::sqrt(rho2);
58 
59  Matrix45 jacobian;
60  jacobian(0, 1) = pred[0] * pred[2];
61  jacobian(0, 4) = rho;
62  jacobian(1, 3) = 1.;
63  jacobian(2, 1) = rho2;
64  jacobian(3, 2) = 1.;
65 
66  return Similarity(jacobian, error.matrix());
67 }
68 
70  const GlobalPoint &priorPosition,
71  const GlobalError &priorError,
72  const GlobalVector &direction,
73  const GlobalError &directionError)
74 {
75  using namespace ROOT::Math;
76 
77  double perp2 = direction.perp2();
78  GlobalVector dir = direction / std::sqrt(perp2);
79  double tip = priorPosition.y() * dir.x() - priorPosition.x() * dir.y();
80  double l = priorPosition.x() * dir.x() + priorPosition.y() * dir.y();
81 
82  prediction_[0] = priorPosition.z() - dir.z() * l;
83  prediction_[1] = tip;
84  prediction_[2] = dir.z();
85  prediction_[3] = dir.phi();
86 
87  Matrix46 jacobian;
88  jacobian(0, 0) = -dir.x() * dir.z();
89  jacobian(1, 0) = -dir.y();
90  jacobian(0, 1) = -dir.y() * dir.z();
91  jacobian(1, 1) = dir.x();
92  jacobian(0, 2) = 1.;
93  jacobian(0, 3) = -dir.z() * priorPosition.x();
94  jacobian(1, 3) = priorPosition.y();
95  jacobian(3, 3) = -dir.y();
96  jacobian(0, 4) = -dir.z() * priorPosition.y();
97  jacobian(1, 4) = -priorPosition.x();
98  jacobian(3, 4) = dir.x();
99  jacobian(0, 5) = -l;
100  jacobian(2, 5) = 1.;
101 
102  Matrix6S origCov;
103  origCov.Place_at(priorError.matrix_new(), 0, 0);
104  origCov.Place_at(directionError.matrix_new() / perp2, 3, 3);
105 
106  covariance_ = Similarity(jacobian, origCov);
107 }
108 
110  const GlobalPoint &priorPosition,
111  const GlobalError &priorError,
112  const GlobalVector &direction,
113  double coneRadius)
114 {
115  double dTheta = std::cosh((double)direction.eta()) * coneRadius;
116 
117  double r2 = direction.mag2();
118  double r = std::sqrt(r2);
119  double perp = direction.perp();
120  double P = direction.x() / perp;
121  double p = direction.y() / perp;
122  double T = direction.z() / r;
123  double t = perp / r;
124  double h2 = dTheta * dTheta;
125  double d2 = coneRadius * coneRadius;
126 
127  GlobalError cov(r2 * (T*T * P*P * h2 + t*t * p*p * d2),
128  r2 * p*P * (T*T * h2 - t*t * d2),
129  r2 * (T*T * p*p * h2 + t*t * P*P * d2),
130  -r2 * t*T * P * h2,
131  -r2 * t*T * p * h2,
132  r2 * t*t * h2);
133 
134  init(priorPosition, priorError, direction, cov);
135 }
136 
138  const CurvilinearTrajectoryParameters &trajectory,
140  prediction_(convert(trajectory)),
142 {
143 }
144 
146  const GlobalTrajectoryParameters &trajectory,
148  prediction_(convert(trajectory)),
149  covariance_(convert(prediction_, error))
150 {
151 }
152 
154  prediction_(convert(
156  GlobalPoint(track.vx(), track.vy(), track.vz()),
157  GlobalVector(track.px(), track.py(), track.pz()),
158  0, 0))),
159  covariance_(convert(prediction_, track.covariance()))
160 {
161 }
162 
164 {
165  using namespace ROOT::Math;
166 
167  double x = std::cos(phi());
168  double y = std::sin(phi());
169 
170  Matrix34 jacobian;
171  jacobian(0, 1) = -y;
172  jacobian(0, 3) = -y * lambda - x * ip();
173  jacobian(1, 1) = x;
174  jacobian(1, 3) = x * lambda - y * ip();
175  jacobian(2, 0) = 1.;
176  jacobian(2, 2) = lambda;
177 
178  return Similarity(jacobian, covariance());
179 }
180 
182 {
183  using namespace ROOT::Math;
184 
185  double x = std::cos(phi());
186  double y = std::sin(phi());
187 
188  Matrix64 jacobian;
189  jacobian(0, 1) = -y;
190  jacobian(0, 3) = -y * lambda - x * ip();
191  jacobian(1, 1) = x;
192  jacobian(1, 3) = x * lambda - y * ip();
193  jacobian(2, 0) = 1.;
194  jacobian(2, 2) = lambda;
195  jacobian(3, 3) = -y;
196  jacobian(4, 3) = x;
197  jacobian(5, 2) = 1.;
198 
199  return Similarity(jacobian, covariance());
200 }
201 
204 {
205  return CurvilinearTrajectoryParameters(0., std::atan(cotTheta()),
206  phi(), ip(), sz(), false);
207 }
208 
210  const MagneticField *fieldProvider) const
211 {
213  0, fieldProvider);
214 }
215 
217 {
218  double rho2I = 1. / rho2();
219  double rhoI = std::sqrt(rho2I);
220 
221  Matrix54 jacobian;
222  jacobian(1, 2) = rho2I;
223  jacobian(2, 3) = 1.;
224  jacobian(3, 1) = 1.;
225  jacobian(4, 0) = rhoI;
226  jacobian(4, 2) = - z() * rhoI * cotTheta() * rho2I;
227 
228  AlgebraicSymMatrix55 err = Similarity(jacobian, covariance());
229  err(0, 0) = 1.;
230 
231  return CurvilinearTrajectoryError(err);
232 }
233 
235  const MagneticField *fieldProvider) const
236 {
237  return FreeTrajectoryState(globalTrajectory(fieldProvider),
238  curvilinearError());
239 }
240 
241 Track GhostTrackPrediction::track(double ndof, double chi2) const
242 {
243  GlobalPoint origin = this->origin();
244  GlobalVector dir = direction().unit();
245 
246  Track::Point point(origin.x(), origin.y(), origin.z());
247  Track::Vector vector(dir.x(), dir.y(), dir.z());
248 
249  return Track(chi2, ndof, point, vector, 0, curvilinearError());
250 }
const GlobalPoint origin() const
ROOT::Math::SMatrix< double, 4, 4, ROOT::Math::MatRepSym< double, 4 > > Error
double lambda(const GlobalPoint &point) const
T perp2() const
Squared magnitude of transverse component.
T perp() const
Magnitude of transverse component.
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Definition: DDAxes.h:10
FreeTrajectoryState fts(const MagneticField *fieldProvider) const
#define P
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
A arg
Definition: Factorize.h:36
void convert(uint32 i, char_uint32 v)
Definition: MsgTools.h:46
const AlgebraicSymMatrix33 & matrix_new() const
GlobalTrajectoryParameters globalTrajectory(const MagneticField *fieldProvider) const
T sqrt(T t)
Definition: SSEVec.h:28
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
void init(const GlobalPoint &priorPosition, const GlobalError &priorError, const GlobalVector &direction, const GlobalError &directionError)
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
math::XYZPoint Point
point in the space
Definition: TrackBase.h:75
CurvilinearTrajectoryError curvilinearError() const
CartesianError cartesianError(double lambda=0.) const
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< float >, ROOT::Math::GlobalCoordinateSystemTag > GlobalVector
vector in glovbal coordinate system
Definition: Vector3D.h:28
Track track(double ndof=0., double chi2=0.) const
CurvilinearTrajectoryParameters curvilinearTrajectory() const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float >, ROOT::Math::GlobalCoordinateSystemTag > GlobalPoint
point in global coordinate system
Definition: Point3D.h:18
const Error & covariance() const
const AlgebraicSymMatrix55 & matrix() const
Square< F >::type sqr(const F &f)
Definition: Square.h:13
dbl *** dir
Definition: mlp_gen.cc:35
ROOT::Math::SVector< double, 4 > Vector
math::XYZVector Vector
spatial vector
Definition: TrackBase.h:73
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
GlobalError positionError(double lambda=0.) const
const GlobalVector direction() const