Fork me on GitHub

source: git/external/TrackCovariance/TrkUtil.cc@ d612dec

Last change on this file since d612dec was ebf40fd, checked in by Franco BEDESCHI <bed@…>, 3 years ago

Major update to handle highly displaced tracks
  • Property mode set to 100644
File size: 23.3 KB
RevLine 
[82db145]1#include "TrkUtil.h"
2#include <iostream>
[c5696dd]3#include <algorithm>
[65776c0]4#include <TSpline.h>
[ebf40fd]5#include <TDecompChol.h>
[82db145]6
7// Constructor
8TrkUtil::TrkUtil(Double_t Bz)
9{
10 fBz = Bz;
11 fGasSel = 0; // Default is He-Isobuthane (90-10)
12 fRmin = 0.0; // Lower DCH radius
13 fRmax = 0.0; // Higher DCH radius
14 fZmin = 0.0; // Lower DCH z
15 fZmax = 0.0; // Higher DCH z
16}
17TrkUtil::TrkUtil()
18{
19 fBz = 0.0;
20 fGasSel = 0; // Default is He-Isobuthane (90-10)
21 fRmin = 0.0; // Lower DCH radius
22 fRmax = 0.0; // Higher DCH radius
23 fZmin = 0.0; // Lower DCH z
24 fZmax = 0.0; // Higher DCH z
25}
26//
27// Destructor
28TrkUtil::~TrkUtil()
29{
30 fBz = 0.0;
31 fGasSel = 0; // Default is He-Isobuthane (90-10)
32 fRmin = 0.0; // Lower DCH radius
33 fRmax = 0.0; // Higher DCH radius
34 fZmin = 0.0; // Lower DCH z
35 fZmax = 0.0; // Higher DCH z
36}
37//
[ebf40fd]38// Distance between two lines
39//
40void TrkUtil::LineDistance(TVector3 x0, TVector3 y0, TVector3 dirx, TVector3 diry, Double_t &sx, Double_t &sy, Double_t &distance)
41{
42 TMatrixDSym M(2);
43 M(0,0) = dirx.Mag2();
44 M(1,1) = diry.Mag2();
45 M(0,1) = -dirx.Dot(diry);
46 M(1,0) = M(0,1);
47 M.Invert();
48 TVectorD c(2);
49 c(0) = dirx.Dot(y0-x0);
50 c(1) = diry.Dot(x0-y0);
51 TVectorD st = M*c;
52 //
53 // Fill output
54 sx = st(0);
55 sy = st(1);
56 //
57 TVector3 x = x0+sx*dirx;
58 TVector3 y = y0+sy*diry;
59 TVector3 d = x-y;
60 distance = d.Mag();
61}
62//
63// Covariance smearing
64//
65TVectorD TrkUtil::CovSmear(TVectorD x, TMatrixDSym C)
66{
67 //
68 // Check arrays
69 //
70 // Consistency of dimensions
71 Int_t Nvec = x.GetNrows();
72 Int_t Nmat = C.GetNrows();
73 if (Nvec != Nmat || Nvec == 0)
74 {
75 std::cout << "TrkUtil::CovSmear: vector/matrix mismatch. Aborting." << std::endl;
76 exit(EXIT_FAILURE);
77 }
78 // Positive diagonal elements
79 for (Int_t i = 0; i < Nvec; i++)
80 {
81 if (C(i, i) <= 0.0)
82 {
83 std::cout << "TrkUtil::CovSmear: covariance matrix has negative diagonal elements. Aborting." << std::endl;
84 exit(EXIT_FAILURE);
85 }
86 }
87 //
88 // Do a Choleski decomposition and random number extraction, with appropriate stabilization
89 //
90 TMatrixDSym CvN = C;
91 TMatrixDSym DCv(Nvec); DCv.Zero();
92 TMatrixDSym DCvInv(Nvec); DCvInv.Zero();
93 for (Int_t id = 0; id < Nvec; id++)
94 {
95 Double_t dVal = TMath::Sqrt(C(id, id));
96 DCv(id, id) = dVal;
97 DCvInv(id, id) = 1.0 / dVal;
98 }
99 CvN.Similarity(DCvInv); // Normalize diagonal to 1
100 TDecompChol Chl(CvN);
101 Bool_t OK = Chl.Decompose(); // Choleski decomposition of normalized matrix
102 if (!OK)
103 {
104 std::cout << "TrkUtil::CovSmear: covariance matrix is not positive definite. Aborting." << std::endl;
105 exit(EXIT_FAILURE);
106 }
107 TMatrixD U = Chl.GetU(); // Get Upper triangular matrix
108 TMatrixD Ut(TMatrixD::kTransposed, U); // Transposed of U (lower triangular)
109 TVectorD r(Nvec);
110 for (Int_t i = 0; i < Nvec; i++)r(i) = gRandom->Gaus(0.0, 1.0); // Array of normal random numbers
111 TVectorD xOut = x + DCv * (Ut * r); // Observed parameter vector
112 //
113 return xOut;
114}
115//
[82db145]116// Helix parameters from position and momentum
117// static
118TVectorD TrkUtil::XPtoPar(TVector3 x, TVector3 p, Double_t Q, Double_t Bz)
119{
120 //
121 TVectorD Par(5);
122 // Transverse parameters
123 Double_t a = -Q * Bz * cSpeed(); // Units are Tesla, GeV and meters
124 Double_t pt = p.Pt();
125 Double_t C = a / (2 * pt); // Half curvature
126 //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl;
[65776c0]127 Double_t r2 = x(0) * x(0) + x(1) * x(1);
[82db145]128 Double_t cross = x(0) * p(1) - x(1) * p(0);
[ebf40fd]129 Double_t T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2);
130 Double_t phi0 = TMath::ATan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); // Phi0
[82db145]131 Double_t D; // Impact parameter D
132 if (pt < 10.0) D = (T - pt) / a;
133 else D = (-2 * cross + a * r2) / (T + pt);
134 //
135 Par(0) = D; // Store D
136 Par(1) = phi0; // Store phi0
137 Par(2) = C; // Store C
138 //Longitudinal parameters
[ebf40fd]139 Double_t B = C * TMath::Sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));
140 Double_t st = TMath::ASin(B) / C;
[82db145]141 Double_t ct = p(2) / pt;
[65776c0]142 Double_t z0;
143 Double_t dot = x(0) * p(0) + x(1) * p(1);
144 if (dot > 0.0) z0 = x(2) - ct * st;
145 else z0 = x(2) + ct * st;
[82db145]146 //
147 Par(3) = z0; // Store z0
148 Par(4) = ct; // Store cot(theta)
149 //
150 return Par;
151}
152// non-static
153TVectorD TrkUtil::XPtoPar(TVector3 x, TVector3 p, Double_t Q)
154{
155 //
156 TVectorD Par(5);
157 Double_t Bz = fBz;
158 Par = XPtoPar(x, p, Q, Bz);
159 //
160 return Par;
161}
162//
163TVector3 TrkUtil::ParToX(TVectorD Par)
164{
165 Double_t D = Par(0);
166 Double_t phi0 = Par(1);
167 Double_t z0 = Par(3);
168 //
169 TVector3 Xval;
[c5696dd]170 Xval(0) = -D * sin(phi0);
171 Xval(1) = D * cos(phi0);
[82db145]172 Xval(2) = z0;
173 //
174 return Xval;
175}
176//
177TVector3 TrkUtil::ParToP(TVectorD Par)
178{
179 if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl;
180 //
[65776c0]181 return ParToP(Par, fBz);
[82db145]182}
183//
184TVector3 TrkUtil::ParToP(TVectorD Par, Double_t Bz)
185{
186 Double_t C = Par(2);
187 Double_t phi0 = Par(1);
188 Double_t ct = Par(4);
189 //
190 TVector3 Pval;
191 Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C);
[c5696dd]192 Pval(0) = pt * cos(phi0);
193 Pval(1) = pt * sin(phi0);
[82db145]194 Pval(2) = pt * ct;
195 //
196 return Pval;
197}
198//
199Double_t TrkUtil::ParToQ(TVectorD Par)
200{
201 return TMath::Sign(1.0, -Par(2));
202}
203
204//
205// Parameter conversion to ACTS format
206TVectorD TrkUtil::ParToACTS(TVectorD Par)
207{
208 TVectorD pACTS(6); // Return vector
209 //
210 Double_t b = -cSpeed() * fBz / 2.;
211 pACTS(0) = 1000 * Par(0); // D from m to mm
212 pACTS(1) = 1000 * Par(3); // z0 from m to mm
213 pACTS(2) = Par(1); // Phi0 is unchanged
[c5696dd]214 pACTS(3) = atan2(1.0, Par(4)); // Theta in [0, pi] range
215 pACTS(4) = Par(2) / (b * sqrt(1 + Par(4) * Par(4))); // q/p in GeV
[82db145]216 pACTS(5) = 0.0; // Time: currently undefined
217 //
218 return pACTS;
219}
220// Covariance conversion to ACTS format
221TMatrixDSym TrkUtil::CovToACTS(TVectorD Par, TMatrixDSym Cov)
222{
223 TMatrixDSym cACTS(6); cACTS.Zero();
224 Double_t b = -cSpeed() * fBz / 2.;
225 //
226 // Fill derivative matrix
227 TMatrixD A(5, 5); A.Zero();
228 Double_t ct = Par(4); // cot(theta)
229 Double_t C = Par(2); // half curvature
230 A(0, 0) = 1000.; // D-D conversion to mm
231 A(1, 2) = 1.0; // phi0-phi0
[c5696dd]232 A(2, 4) = 1.0 / (sqrt(1.0 + ct * ct) * b); // q/p-C
[82db145]233 A(3, 1) = 1000.; // z0-z0 conversion to mm
234 A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta)
235 A(4, 4) = -C * ct / (b * pow(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta)
236 //
237 TMatrixDSym Cv = Cov;
238 TMatrixD At(5, 5);
239 At.Transpose(A);
240 Cv.Similarity(At);
241 TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv;
242 cACTS(5, 5) = 0.1; // Currently undefined: set to arbitrary value to avoid crashes
243 //
244 return cACTS;
245}
246//
247// Parameter conversion to ILC format
248TVectorD TrkUtil::ParToILC(TVectorD Par)
249{
250 TVectorD pILC(5); // Return vector
251 //
252 pILC(0) = Par(0) * 1.0e3; // d0 in mm
253 pILC(1) = Par(1); // phi0 is unchanged
254 pILC(2) = -2 * Par(2) * 1.0e-3; // w in mm^-1
255 pILC(3) = Par(3) * 1.0e3; // z0 in mm
256 pILC(4) = Par(4); // tan(lambda) = cot(theta)
257 //
258 return pILC;
259}
260// Covariance conversion to ILC format
261TMatrixDSym TrkUtil::CovToILC(TMatrixDSym Cov)
262{
263 TMatrixDSym cILC(5); cILC.Zero();
264 //
265 // Fill derivative matrix
266 TMatrixD A(5, 5); A.Zero();
267 //
268 A(0, 0) = 1.0e3; // D-d0 in mm
269 A(1, 1) = 1.0; // phi0-phi0
270 A(2, 2) = -2.0e-3; // w-C
271 A(3, 3) = 1.0e3; // z0-z0 conversion to mm
272 A(4, 4) = 1.0; // tan(lambda) - cot(theta)
273 //
274 TMatrixDSym Cv = Cov;
275 TMatrixD At(5, 5);
276 At.Transpose(A);
277 Cv.Similarity(At);
278 cILC = Cv;
279 //
280 return cILC;
281}
282//
283// Conversion from meters to mm
284TVectorD TrkUtil::ParToMm(TVectorD Par) // Parameter conversion
285{
286 TVectorD Pmm(5); // Return vector
287 //
288 Pmm(0) = Par(0) * 1.0e3; // d0 in mm
289 Pmm(1) = Par(1); // phi0 is unchanged
290 Pmm(2) = Par(2) * 1.0e-3; // C in mm^-1
291 Pmm(3) = Par(3) * 1.0e3; // z0 in mm
292 Pmm(4) = Par(4); // tan(lambda) = cot(theta) unchanged
293 //
294 return Pmm;
295}
296TMatrixDSym TrkUtil::CovToMm(TMatrixDSym Cov) // Covariance conversion
297{
298 TMatrixDSym Cmm(5); Cmm.Zero();
299 //
300 // Fill derivative matrix
301 TMatrixD A(5, 5); A.Zero();
302 //
303 A(0, 0) = 1.0e3; // D-d0 in mm
304 A(1, 1) = 1.0; // phi0-phi0
305 A(2, 2) = 1.0e-3; // C-C
306 A(3, 3) = 1.0e3; // z0-z0 conversion to mm
307 A(4, 4) = 1.0; // lambda - cot(theta)
308 //
309 TMatrixDSym Cv = Cov;
310 TMatrixD At(5, 5);
311 At.Transpose(A);
312 Cv.Similarity(At);
313 Cmm = Cv;
314 //
315 return Cmm;
[ebf40fd]316}//
317// Regularized symmetric matrix inversion
318//
319TMatrixDSym TrkUtil::RegInv(TMatrixDSym& Min)
320{
321 TMatrixDSym M = Min; // Decouple from input
322 Int_t N = M.GetNrows(); // Matrix size
323 TMatrixDSym D(N); D.Zero(); // Normaliztion matrix
324 TMatrixDSym R(N); // Normarized matrix
325 TMatrixDSym Rinv(N); // Inverse of R
326 TMatrixDSym Minv(N); // Inverse of M
327 //
328 // Check for 0's and normalize
329 for (Int_t i = 0; i < N; i++)
330 {
331 if (M(i, i) != 0.0) D(i, i) = 1. / TMath::Sqrt(TMath::Abs(M(i, i)));
332 else D(i, i) = 1.0;
333 }
334 R = M.Similarity(D);
335 //
336 // Recursive algorithms stops when N = 2
337 //
338 //****************
339 // case N = 2 ***
340 //****************
341 if (N == 2)
342 {
343 Double_t det = R(0, 0) * R(1, 1) - R(0, 1) * R(1, 0);
344 if (det == 0)
345 {
346 std::cout << "VertexFit::RegInv: null determinant for N = 2" << std::endl;
347 Rinv.Zero(); // Return null matrix
348 }
349 else
350 {
351 // invert matrix
352 Rinv(0, 0) = R(1, 1);
353 Rinv(0, 1) = -R(0, 1);
354 Rinv(1, 0) = Rinv(0, 1);
355 Rinv(1, 1) = R(0, 0);
356 Rinv *= 1. / det;
357 }
358 }
359 //****************
360 // case N > 2 ***
361 //****************
362 else
363 {
364 // Break up matrix
365 TMatrixDSym Q = R.GetSub(0, N - 2, 0, N - 2); // Upper left
366 TVectorD p(N - 1);
367 for (Int_t i = 0; i < N - 1; i++)p(i) = R(N - 1, i);
368 Double_t q = R(N - 1, N - 1);
369 //Invert pieces and re-assemble
370 TMatrixDSym Ainv(N - 1);
371 TMatrixDSym A(N - 1);
372 if (TMath::Abs(q) > 1.0e-15)
373 {
374 // Case |q| > 0
375 Ainv.Rank1Update(p, -1.0 / q);
376 Ainv += Q;
377 A = RegInv(Ainv); // Recursive call
378 TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A;
379 //
380 TVectorD b = (-1.0 / q) * (A * p);
381 for (Int_t i = 0; i < N - 1; i++)
382 {
383 Rinv(N - 1, i) = b(i);
384 Rinv(i, N - 1) = b(i);
385 }
386 //
387 Double_t pdotb = 0.;
388 for (Int_t i = 0; i < N - 1; i++)pdotb += p(i) * b(i);
389 Double_t c = (1.0 - pdotb) / q;
390 Rinv(N - 1, N - 1) = c;
391 }
392 else
393 {
394 // case q = 0
395 TMatrixDSym Qinv = RegInv(Q); // Recursive call
396 Double_t a = Qinv.Similarity(p);
397 Double_t c = -1.0 / a;
398 Rinv(N - 1, N - 1) = c;
399 //
400 TVectorD b = (1.0 / a) * (Qinv * p);
401 for (Int_t i = 0; i < N - 1; i++)
402 {
403 Rinv(N - 1, i) = b(i);
404 Rinv(i, N - 1) = b(i);
405 }
406 //
407 A.Rank1Update(p, -1 / a);
408 A += Q;
409 A.Similarity(Qinv);
410 TMatrixDSub(Rinv, 0, N - 2, 0, N - 2) = A;
411 }
412 }
413 Minv = Rinv.Similarity(D);
414 return Minv;
415}
416//
417// Track tracjectory
418//
419TVector3 TrkUtil::Xtrack(TVectorD par, Double_t s)
420{
421 //
422 // unpack parameters
423 Double_t D = par(0);
424 Double_t p0 = par(1);
425 Double_t C = par(2);
426 Double_t z0 = par(3);
427 Double_t ct = par(4);
428 //
429 Double_t x = -D * TMath::Sin(p0) + (TMath::Sin(s + p0) - TMath::Sin(p0)) / (2 * C);
430 Double_t y = D * TMath::Cos(p0) - (TMath::Cos(s + p0) - TMath::Cos(p0)) / (2 * C);
431 Double_t z = z0 + ct * s / (2 * C);
432 //
433 TVector3 Xt(x, y, z);
434 return Xt;
435}
436//
437// Track derivatives
438//
439// Constant radius
440// R-Phi
441TVectorD TrkUtil::derRphi_R(TVectorD par, Double_t R)
442{
443 TVectorD dRphi(5); // return vector
444 //
445 // unpack parameters
446 Double_t D = par(0);
447 Double_t C = par(2);
448 //
449 Double_t s = 2 * TMath::ASin(C * TMath::Sqrt((R * R - D * D)/(1 + 2 * C * D)));
450 TVector3 X = Xtrack(par, s); // Intersection point
451 TVector3 v(-X.y()/R, X.x()/R, 0.); // measurement direction
452 TMatrixD derX = derXdPar(par, s); // dX/dp
453 TVectorD derXs = derXds(par, s); // dX/ds
454 TVectorD ders = dsdPar_R(par, R); // ds/dp
455 //
456 for (Int_t i = 0; i < 5; i++)
457 {
458 dRphi(i) = 0.;
459 for (Int_t j = 0; j < 3; j++)
460 {
461 dRphi(i) += v(j) * (derX(j, i) + derXs(j) * ders(i));
462 }
463 }
464 //
465 return dRphi;
466}
467// z
468TVectorD TrkUtil::derZ_R(TVectorD par, Double_t R)
469{
470
471 TVectorD dZ(5); // return vector
472 //
473 // unpack parameters
474 Double_t D = par(0);
475 Double_t C = par(2);
476 //
477 Double_t s = 2 * TMath::ASin(C * TMath::Sqrt((R * R - D * D)/(1 + 2 * C * D))); // phase
478 TVector3 v(0., 0., 1.); // measurement direction
479 TMatrixD derX = derXdPar(par, s); // dX/dp
480 TVectorD derXs = derXds(par, s); // dX/ds
481 TVectorD ders = dsdPar_R(par, R); // ds/dp
482 //
483 for (Int_t i = 0; i < 5; i++)
484 {
485 dZ(i) = 0.;
486 for (Int_t j = 0; j < 3; j++)
487 {
488 dZ(i) += v(j) * (derX(j, i) + derXs(j) * ders(i));
489 }
490 }
491 //
492 return dZ;
493}
494//
495// constant z
496// R-Phi
497TVectorD TrkUtil::derRphi_Z(TVectorD par, Double_t z)
498{
499 TVectorD dRphi(5); // return vector
500 //
501 // unpack parameters
502 Double_t C = par(2);
503 Double_t z0 = par(3);
504 Double_t ct = par(4);
505 //
506 Double_t s = 2 * C * (z - z0) / ct;
507 TVector3 X = Xtrack(par, s); // Intersection point
508 TVector3 v(-X.y() / X.Pt(), X.x() / X.Pt(), 0.); // measurement direction
509 TMatrixD derX = derXdPar(par, s); // dX/dp
510 TVectorD derXs = derXds(par, s); // dX/ds
511 TVectorD ders = dsdPar_z(par, z); // ds/dp
512 //
513 for (Int_t i = 0; i < 5; i++)
514 {
515 dRphi(i) = 0.;
516 for (Int_t j = 0; j < 3; j++)
517 {
518 dRphi(i) += v(j) * (derX(j, i) + derXs(j) * ders(i));
519 }
520 }
521 //
522 return dRphi;
523
524}
525// R
526TVectorD TrkUtil::derR_Z(TVectorD par, Double_t z)
527{
528 TVectorD dR(5); // return vector
529 //
530 // unpack parameters
531 Double_t C = par(2);
532 Double_t z0 = par(3);
533 Double_t ct = par(4);
534 //
535 Double_t s = 2 * C * (z - z0) / ct;
536 TVector3 X = Xtrack(par, s); // Intersection point
537 TVector3 v(X.x() / X.Pt(), X.y() / X.Pt(), 0.); // measurement direction
538 TMatrixD derX = derXdPar(par, s); // dX/dp
539 TVectorD derXs = derXds(par, s); // dX/ds
540 TVectorD ders = dsdPar_z(par, z); // ds/dp
541 //
542 for (Int_t i = 0; i < 5; i++)
543 {
544 dR(i) = 0.;
545 for (Int_t j = 0; j < 3; j++)
546 {
547 dR(i) += v(j) * (derX(j, i) + derXs(j) * ders(i));
548 }
549 }
550 //
551 return dR;
552
553}
554//
555// derivatives of track trajectory
556//
557// dX/dPar
558TMatrixD TrkUtil::derXdPar(TVectorD par, Double_t s)
559{
560 TMatrixD dxdp(3, 5); // return matrix
561 //
562 // unpack parameters
563 Double_t D = par(0);
564 Double_t p0 = par(1);
565 Double_t C = par(2);
566 Double_t z0 = par(3);
567 Double_t ct = par(4);
568 //
569 // derivatives
570 // dx/dD
571 dxdp(0, 0) = -TMath::Sin(p0);
572 dxdp(1, 0) = TMath::Cos(p0);
573 dxdp(2, 0) = 0.;
574 // dx/dphi0
575 dxdp(0, 1) = -D * TMath::Cos(p0) + (TMath::Cos(s + p0) - TMath::Cos(p0)) / (2 * C);
576 dxdp(1, 1) = -D * TMath::Sin(p0) + (TMath::Sin(s + p0) - TMath::Sin(p0)) / (2 * C);
577 dxdp(2, 1) = 0;
578 // dx/dC
579 dxdp(0, 2) = -(TMath::Sin(s + p0) - TMath::Sin(p0)) / (2 * C * C);
580 dxdp(1, 2) = (TMath::Cos(s + p0) - TMath::Cos(p0)) / (2 * C * C);
581 dxdp(2, 2) = -ct * s / (2 * C * C);
582 // dx/dz0
583 dxdp(0, 3) = 0;
584 dxdp(1, 3) = 0;
585 dxdp(2, 3) = 1.;
586 // dx/dCtg
587 dxdp(0, 4) = 0;
588 dxdp(1, 4) = 0;
589 dxdp(2, 4) = s / (2 * C);
590 //
591 return dxdp;
592}
593//
594// dX/ds
595//
596TVectorD TrkUtil::derXds(TVectorD par, Double_t s)
597{
598 TVectorD dxds(3); // return vector
599 //
600 // unpack parameters
601 Double_t p0 = par(1);
602 Double_t C = par(2);
603 Double_t ct = par(4);
604 //
605 // dX/ds
606 dxds(0) = TMath::Cos(s + p0) / (2 * C);
607 dxds(1) = TMath::Sin(s + p0) / (2 * C);
608 dxds(2) = ct / (2 * C);
609 //
610 return dxds;
611}
612//
613// derivative of trajectory phase s
614//Constant R
615TVectorD TrkUtil::dsdPar_R(TVectorD par, Double_t R)
616{
617 TVectorD dsdp(5); // return vector
618 //
619 // unpack parameters
620 Double_t D = par(0);
621 Double_t p0 = par(1);
622 Double_t C = par(2);
623 //
624 // derivatives
625 Double_t opCD = 1. + 2 * C * D;
626 Double_t A = C*TMath::Sqrt((R*R-D*D)/opCD);
627 Double_t sqA0 = TMath::Sqrt(1. - A * A);
628 Double_t dMin = 0.01;
629 Double_t sqA = TMath::Max(dMin, sqA0); // Protect against divergence
630 //
631 dsdp(0) = -2 * C * C * (D * (1. + C * D) + C * R * R) / (A * sqA * opCD * opCD);
632 dsdp(1) = 0;
633 dsdp(2) = 2 * A * (1 + C * D) / (C * sqA * opCD);
634 dsdp(3) = 0;
635 dsdp(4) = 0;
636 //
637 return dsdp;
638}
639// Constant z
640TVectorD TrkUtil::dsdPar_z(TVectorD par, Double_t z)
641{
642 TVectorD dsdp(5); // return vector
643 //
644 // unpack parameters
645 Double_t C = par(2);
646 Double_t z0 = par(3);
647 Double_t ct = par(4);
648 //
649 // derivatives
650 //
651 dsdp(0) = 0;
652 dsdp(1) = 0;
653 dsdp(2) = 2*(z-z0)/ct;
654 dsdp(3) = -2*C/ct;
655 dsdp(4) = -2*C*(z-z0)/(ct*ct);
656 //
657 return dsdp;
[82db145]658}
659//
660// Setup chamber volume
661void TrkUtil::SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax)
662{
663 fRmin = Rmin; // Lower DCH radius
664 fRmax = Rmax; // Higher DCH radius
665 fZmin = Zmin; // Lower DCH z
666 fZmax = Zmax; // Higher DCH z
667}
668//
669// Get Trakck length inside DCH volume
670Double_t TrkUtil::TrkLen(TVectorD Par)
671{
672 Double_t tLength = 0.0;
673 // Check if geometry is initialized
674 if (fZmin == 0.0 && fZmax == 0.0)
675 {
676 // No geometry set so send a warning and return 0
677 std::cout << "TrkUtil::TrkLen() called without a DCH volume defined" << std::endl;
678 }
679 else
680 {
681 //******************************************************************
682 // Determine the track length inside the chamber ****
683 //******************************************************************
684 //
685 // Track pararameters
686 Double_t D = Par(0); // Transverse impact parameter
687 Double_t phi0 = Par(1); // Transverse direction at minimum approach
688 Double_t C = Par(2); // Half curvature
689 Double_t z0 = Par(3); // Z at minimum approach
690 Double_t ct = Par(4); // cot(theta)
691 //std::cout << "TrkUtil:: parameters: D= " << D << ", phi0= " << phi0
692 // << ", C= " << C << ", z0= " << z0 << ", ct= " << ct << std::endl;
693 //
[65776c0]694 // Track length per unit phase change
695 Double_t Scale = sqrt(1.0 + ct * ct) / (2.0 * TMath::Abs(C));
[82db145]696 //
697 // Find intersections with chamber boundaries
698 //
[65776c0]699 Double_t phRin = 0.0; // phase of inner cylinder
700 Double_t phRin2 = 0.0; // phase of inner cylinder intersection (2nd branch)
[82db145]701 Double_t phRhi = 0.0; // phase of outer cylinder intersection
702 Double_t phZmn = 0.0; // phase of left wall intersection
703 Double_t phZmx = 0.0; // phase of right wall intersection
704 // ... with inner cylinder
[65776c0]705 Double_t Rtop = TMath::Abs((1.0 + C * D) / C);
[82db145]706
707 if (Rtop > fRmin && TMath::Abs(D) < fRmin) // *** don't treat large D tracks for the moment ***
708 {
[65776c0]709 Double_t ph = 2 * asin(C * sqrt((fRmin * fRmin - D * D) / (1.0 + 2.0 * C * D)));
710 Double_t z = z0 + ct * ph / (2.0 * C);
[82db145]711
712 //std::cout << "Rin intersection: ph = " << ph<<", z= "<<z << std::endl;
713
[65776c0]714 if (z < fZmax && z > fZmin) phRin = TMath::Abs(ph); // Intersection inside chamber volume
[82db145]715 //
716 // Include second branch of loopers
[59ba063]717 Double_t Pi = 3.14159265358979323846;
[65776c0]718 Double_t ph2 = 2 * Pi - TMath::Abs(ph);
[82db145]719 if (ph < 0)ph2 = -ph2;
720 z = z0 + ct * ph2 / (2.0 * C);
721 if (z < fZmax && z > fZmin) phRin2 = TMath::Abs(ph2); // Intersection inside chamber volume
722 }
723 // ... with outer cylinder
724 if (Rtop > fRmax && TMath::Abs(D) < fRmax) // *** don't treat large D tracks for the moment ***
725 {
[65776c0]726 Double_t ph = 2 * asin(C * sqrt((fRmax * fRmax - D * D) / (1.0 + 2.0 * C * D)));
727 Double_t z = z0 + ct * ph / (2.0 * C);
728 if (z < fZmax && z > fZmin) phRhi = TMath::Abs(ph); // Intersection inside chamber volume
[82db145]729 }
730 // ... with left wall
731 Double_t Zdir = (fZmin - z0) / ct;
732 if (Zdir > 0.0)
733 {
[65776c0]734 Double_t ph = 2.0 * C * Zdir;
735 Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C));
736 if (Rint < fRmax && Rint > fRmin) phZmn = TMath::Abs(ph); // Intersection inside chamber volume
[82db145]737 }
738 // ... with right wall
739 Zdir = (fZmax - z0) / ct;
740 if (Zdir > 0.0)
741 {
[65776c0]742 Double_t ph = 2.0 * C * Zdir;
743 Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C));
744 if (Rint < fRmax && Rint > fRmin) phZmx = TMath::Abs(ph); // Intersection inside chamber volume
[82db145]745 }
746 //
747 // Order phases and keep the lowest two non-zero ones
748 //
749 const Int_t Nint = 5;
750 Double_t dPhase = 0.0; // Phase difference between two close intersections
751 Double_t ph_arr[Nint] = { phRin, phRin2, phRhi, phZmn, phZmx };
[59ba063]752 std::sort(ph_arr, ph_arr + Nint);
[82db145]753 Int_t iPos = -1; // First element > 0
754 for (Int_t i = 0; i < Nint; i++)
755 {
[59ba063]756 if (ph_arr[i] <= 0.0) iPos = i;
[82db145]757 }
758
759 if (iPos < Nint - 2)
760 {
[59ba063]761 dPhase = ph_arr[iPos + 2] - ph_arr[iPos + 1];
[65776c0]762 tLength = dPhase * Scale;
[82db145]763 }
764 }
765 return tLength;
766}
767//
768// Return number of ionization clusters
[65776c0]769Bool_t TrkUtil::IonClusters(Double_t& Ncl, Double_t mass, TVectorD Par)
[82db145]770{
771 //
772 // Units are meters/Tesla/GeV
773 //
774 Ncl = 0.0;
775 Bool_t Signal = kFALSE;
776 Double_t tLen = 0;
777 // Check if geometry is initialized
778 if (fZmin == 0.0 && fZmax == 0.0)
779 {
780 // No geometry set so send a warning and return 0
781 std::cout << "TrkUtil::IonClusters() called without a volume defined" << std::endl;
782 }
783 else tLen = TrkLen(Par);
784
785 //******************************************************************
786 // Now get the number of clusters ****
787 //******************************************************************
788 //
789 Double_t muClu = 0.0; // mean number of clusters
790 Double_t bg = 0.0; // beta*gamma
791 Ncl = 0.0;
792 if (tLen > 0.0)
793 {
794 Signal = kTRUE;
795 //
796 // Find beta*gamma
797 if (fBz == 0.0)
798 {
799 Signal = kFALSE;
800 std::cout << "TrkUtil::IonClusters: Please set Bz!!!" << std::endl;
801 }
802 else
803 {
804 TVector3 p = ParToP(Par);
805 bg = p.Mag() / mass;
[65776c0]806 muClu = Nclusters(bg) * tLen; // Avg. number of clusters
[82db145]807
808 Ncl = gRandom->PoissonD(muClu); // Actual number of clusters
809 }
810
811 }
[65776c0]812 //
[82db145]813 return Signal;
814}
815//
816//
[a95da74]817Double_t TrkUtil::Nclusters(Double_t begam)
[82db145]818{
819 Int_t Opt = fGasSel;
820 Double_t Nclu = Nclusters(begam, Opt);
821 //
822 return Nclu;
823}
824//
825Double_t TrkUtil::Nclusters(Double_t begam, Int_t Opt) {
826 //
827 // Opt = 0: He 90 - Isobutane 10
828 // = 1: pure He
829 // = 2: Argon 50 - Ethane 50
830 // = 3: pure Argon
831 //
832 //
[4df491e]833 const Int_t Npt = 18;
834 Double_t bg[Npt] = { 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10.,
835 12., 15., 20., 50., 100., 200., 500., 1000., 10000. };
836 //
837 // He 90 - Isobutane 10
838 Double_t ncl_He_Iso[Npt] = { 42.94, 23.6,18.97,12.98,12.2,12.13,
839 12.24,12.73,13.03,13.29,13.63,14.08,15.56,16.43,16.8,16.95,16.98, 16.98 };
840 //
841 // pure He
842 Double_t ncl_He[Npt] = { 11.79,6.5,5.23,3.59,3.38,3.37,3.4,3.54,3.63,
843 3.7,3.8,3.92,4.33,4.61,4.78,4.87,4.89, 4.89 };
844 //
845 // Argon 50 - Ethane 50
846 Double_t ncl_Ar_Eth[Npt] = { 130.04,71.55,57.56,39.44,37.08,36.9,
847 37.25,38.76,39.68,40.49,41.53,42.91,46.8,48.09,48.59,48.85,48.93,48.93 };
848 //
849 // pure Argon
850 Double_t ncl_Ar[Npt] = { 88.69,48.93,39.41,27.09,25.51,25.43,25.69,
851 26.78,27.44,28.02,28.77,29.78,32.67,33.75,34.24,34.57,34.68, 34.68 };
852 //
853 Double_t ncl[Npt];
[65776c0]854 switch (Opt)
855 {
856 case 0: std::copy(ncl_He_Iso, ncl_He_Iso + Npt, ncl); // He-Isobutane
[a95da74]857 break;
[65776c0]858 case 1: std::copy(ncl_He, ncl_He + Npt, ncl); // pure He
[4df491e]859 break;
[65776c0]860 case 2: std::copy(ncl_Ar_Eth, ncl_Ar_Eth + Npt, ncl); // Argon - Ethane
[4df491e]861 break;
[65776c0]862 case 3: std::copy(ncl_Ar, ncl_Ar + Npt, ncl); // pure Argon
[4df491e]863 break;
[65776c0]864 }
865 //
866 Double_t interp = 0.0;
867 TSpline3* sp3 = new TSpline3("sp3", bg, ncl, Npt);
868 if (begam > bg[0] && begam < bg[Npt - 1]) interp = sp3->Eval(begam);
869 return 100 * interp;
[82db145]870}
871//
[65776c0]872Double_t TrkUtil::funcNcl(Double_t* xp, Double_t* par) {
[82db145]873 Double_t bg = xp[0];
874 return Nclusters(bg);
875}
876//
877void TrkUtil::SetGasMix(Int_t Opt)
878{
879 if (Opt < 0 || Opt > 3)
880 {
881 std::cout << "TrkUtil::SetGasMix Gas option not allowed. No action."
882 << std::endl;
883 }
884 else fGasSel = Opt;
[c5696dd]885}
Note: See TracBrowser for help on using the repository browser.