Fork me on GitHub

source: git/external/TrackCovariance/TrkUtil.cc@ 3dc190c

Last change on this file since 3dc190c was 92b8d11, checked in by Franco BEDESCHI <bed@…>, 4 years ago

Fix Sqrt for ROOT5
  • Property mode set to 100644
File size: 5.1 KB
Line 
1#include "TrkUtil.h"
2#include <iostream>
3
4// Constructor
5TrkUtil::TrkUtil(Double_t Bz)
6{
7 fBz = Bz;
8}
9TrkUtil::TrkUtil()
10{
11 fBz = 0.0;
12}
13//
14// Destructor
15TrkUtil::~TrkUtil()
16{
17 fBz = 0.0;
18}
19//
20// Helix parameters from position and momentum
21// static
22TVectorD TrkUtil::XPtoPar(TVector3 x, TVector3 p, Double_t Q, Double_t Bz)
23{
24 //
25 TVectorD Par(5);
26 // Transverse parameters
27 Double_t a = -Q * Bz * cSpeed(); // Units are Tesla, GeV and meters
28 Double_t pt = p.Pt();
29 Double_t C = a / (2 * pt); // Half curvature
30 //cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << endl;
31 Double_t r2 = x.Perp2();
32 Double_t cross = x(0) * p(1) - x(1) * p(0);
33 Double_t T = sqrt(pt * pt - 2 * a * cross + a * a * r2);
34 Double_t phi0 = atan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); // Phi0
35 Double_t D; // Impact parameter D
36 if (pt < 10.0) D = (T - pt) / a;
37 else D = (-2 * cross + a * r2) / (T + pt);
38 //
39 Par(0) = D; // Store D
40 Par(1) = phi0; // Store phi0
41 Par(2) = C; // Store C
42 //Longitudinal parameters
43 Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));
44 Double_t st = asin(B) / C;
45 Double_t ct = p(2) / pt;
46 Double_t z0 = x(2) - ct * st;
47 //
48 Par(3) = z0; // Store z0
49 Par(4) = ct; // Store cot(theta)
50 //
51 return Par;
52}
53// non-static
54TVectorD TrkUtil::XPtoPar(TVector3 x, TVector3 p, Double_t Q)
55{
56 //
57 TVectorD Par(5);
58 Double_t Bz = fBz;
59 Par = XPtoPar(x, p, Q, Bz);
60 //
61 return Par;
62}
63//
64TVector3 TrkUtil::ParToX(TVectorD Par)
65{
66 Double_t D = Par(0);
67 Double_t phi0 = Par(1);
68 Double_t z0 = Par(3);
69 //
70 TVector3 Xval;
71 Xval(0) = -D * sin(phi0);
72 Xval(1) = D * cos(phi0);
73 Xval(2) = z0;
74 //
75 return Xval;
76}
77//
78TVector3 TrkUtil::ParToP(TVectorD Par)
79{
80 if (fBz == 0.0)
81std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl;
82 //
83 return ParToP(Par,fBz);
84}
85//
86TVector3 TrkUtil::ParToP(TVectorD Par, Double_t Bz)
87{
88 Double_t C = Par(2);
89 Double_t phi0 = Par(1);
90 Double_t ct = Par(4);
91 //
92 TVector3 Pval;
93 Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C);
94 Pval(0) = pt * cos(phi0);
95 Pval(1) = pt * sin(phi0);
96 Pval(2) = pt * ct;
97 //
98 return Pval;
99}
100//
101Double_t TrkUtil::ParToQ(TVectorD Par)
102{
103 return TMath::Sign(1.0, -Par(2));
104}
105
106//
107// Parameter conversion to ACTS format
108TVectorD TrkUtil::ParToACTS(TVectorD Par)
109{
110 TVectorD pACTS(6); // Return vector
111 //
112 Double_t b = -cSpeed() * fBz / 2.;
113 pACTS(0) = 1000 * Par(0); // D from m to mm
114 pACTS(1) = 1000 * Par(3); // z0 from m to mm
115 pACTS(2) = Par(1); // Phi0 is unchanged
116 pACTS(3) = atan2(1.0, Par(4)); // Theta in [0, pi] range
117 pACTS(4) = Par(2) / (b * sqrt(1 + Par(4) * Par(4))); // q/p in GeV
118 pACTS(5) = 0.0; // Time: currently undefined
119 //
120 return pACTS;
121}
122// Covariance conversion to ACTS format
123TMatrixDSym TrkUtil::CovToACTS(TVectorD Par, TMatrixDSym Cov)
124{
125 TMatrixDSym cACTS(6); cACTS.Zero();
126 Double_t b = -cSpeed() * fBz / 2.;
127 //
128 // Fill derivative matrix
129 TMatrixD A(5, 5); A.Zero();
130 Double_t ct = Par(4); // cot(theta)
131 Double_t C = Par(2); // half curvature
132 A(0, 0) = 1000.; // D-D conversion to mm
133 A(1, 2) = 1.0; // phi0-phi0
134 A(2, 4) = 1.0 / (sqrt(1.0 + ct * ct) * b); // q/p-C
135 A(3, 1) = 1000.; // z0-z0 conversion to mm
136 A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta)
137 A(4, 4) = -C * ct / (b * pow(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta)
138 //
139 TMatrixDSym Cv = Cov;
140 TMatrixD At(5, 5);
141 At.Transpose(A);
142 Cv.Similarity(At);
143 TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv;
144 cACTS(5, 5) = 0.1; // Currently undefined: set to arbitrary value to avoid crashes
145 //
146 return cACTS;
147}
148//
149// Parameter conversion to ILC format
150TVectorD TrkUtil::ParToILC(TVectorD Par)
151{
152 TVectorD pILC(5); // Return vector
153 //
154 pILC(0) = Par(0) * 1.0e3; // d0 in mm
155 pILC(1) = Par(1); // phi0 is unchanged
156 pILC(2) = -2 * Par(2) * 1.0e-3; // w in mm^-1
157 pILC(3) = Par(3) * 1.0e3; // z0 in mm
158 pILC(4) = Par(4); // tan(lambda) = cot(theta)
159 //
160 return pILC;
161}
162// Covariance conversion to ILC format
163TMatrixDSym TrkUtil::CovToILC(TMatrixDSym Cov)
164{
165 TMatrixDSym cILC(5); cILC.Zero();
166 //
167 // Fill derivative matrix
168 TMatrixD A(5, 5); A.Zero();
169 //
170 A(0, 0) = 1.0e3; // D-d0 in mm
171 A(1, 1) = 1.0; // phi0-phi0
172 A(2, 2) = -2.0e-3; // w-C
173 A(3, 3) = 1.0e3; // z0-z0 conversion to mm
174 A(4, 4) = 1.0; // tan(lambda) - cot(theta)
175 //
176 TMatrixDSym Cv = Cov;
177 TMatrixD At(5, 5);
178 At.Transpose(A);
179 Cv.Similarity(At);
180 cILC = Cv;
181 //
182 return cILC;
183}
184//
185// Conversion from meters to mm
186TVectorD TrkUtil::ParToMm(TVectorD Par) // Parameter conversion
187{
188 TVectorD Pmm(5); // Return vector
189 //
190 Pmm(0) = Par(0) * 1.0e3; // d0 in mm
191 Pmm(1) = Par(1); // phi0 is unchanged
192 Pmm(2) = Par(2) * 1.0e-3; // C in mm^-1
193 Pmm(3) = Par(3) * 1.0e3; // z0 in mm
194 Pmm(4) = Par(4); // tan(lambda) = cot(theta) unchanged
195 //
196 return Pmm;
197}
198TMatrixDSym TrkUtil::CovToMm(TMatrixDSym Cov) // Covariance conversion
199{
200 TMatrixDSym Cmm(5); Cmm.Zero();
201 //
202 // Fill derivative matrix
203 TMatrixD A(5, 5); A.Zero();
204 //
205 A(0, 0) = 1.0e3; // D-d0 in mm
206 A(1, 1) = 1.0; // phi0-phi0
207 A(2, 2) = 1.0e-3; // C-C
208 A(3, 3) = 1.0e3; // z0-z0 conversion to mm
209 A(4, 4) = 1.0; // lambda - cot(theta)
210 //
211 TMatrixDSym Cv = Cov;
212 TMatrixD At(5, 5);
213 At.Transpose(A);
214 Cv.Similarity(At);
215 Cmm = Cv;
216 //
217 return Cmm;
218}
Note: See TracBrowser for help on using the repository browser.