Fork me on GitHub

source: git/external/TrackCovariance/TrkUtil.cc@ 27197df

Last change on this file since 27197df was 4df491e, checked in by Franco BEDESCHI <bed@…>, 4 years ago

Fixed cluster counting inits
  • Property mode set to 100644
File size: 15.1 KB
RevLine 
[82db145]1#include "TrkUtil.h"
2#include <iostream>
[c5696dd]3#include <algorithm>
[82db145]4
5// Constructor
6TrkUtil::TrkUtil(Double_t Bz)
7{
8 fBz = Bz;
9 fGasSel = 0; // Default is He-Isobuthane (90-10)
10 fRmin = 0.0; // Lower DCH radius
11 fRmax = 0.0; // Higher DCH radius
12 fZmin = 0.0; // Lower DCH z
13 fZmax = 0.0; // Higher DCH z
14}
15TrkUtil::TrkUtil()
16{
17 fBz = 0.0;
18 fGasSel = 0; // Default is He-Isobuthane (90-10)
19 fRmin = 0.0; // Lower DCH radius
20 fRmax = 0.0; // Higher DCH radius
21 fZmin = 0.0; // Lower DCH z
22 fZmax = 0.0; // Higher DCH z
23}
24//
25// Destructor
26TrkUtil::~TrkUtil()
27{
28 fBz = 0.0;
29 fGasSel = 0; // Default is He-Isobuthane (90-10)
30 fRmin = 0.0; // Lower DCH radius
31 fRmax = 0.0; // Higher DCH radius
32 fZmin = 0.0; // Lower DCH z
33 fZmax = 0.0; // Higher DCH z
34}
35//
36// Helix parameters from position and momentum
37// static
38TVectorD TrkUtil::XPtoPar(TVector3 x, TVector3 p, Double_t Q, Double_t Bz)
39{
40 //
41 TVectorD Par(5);
42 // Transverse parameters
43 Double_t a = -Q * Bz * cSpeed(); // Units are Tesla, GeV and meters
44 Double_t pt = p.Pt();
45 Double_t C = a / (2 * pt); // Half curvature
46 //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl;
47 Double_t r2 = x.Perp2();
48 Double_t cross = x(0) * p(1) - x(1) * p(0);
[c5696dd]49 Double_t T = sqrt(pt * pt - 2 * a * cross + a * a * r2);
50 Double_t phi0 = atan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); // Phi0
[82db145]51 Double_t D; // Impact parameter D
52 if (pt < 10.0) D = (T - pt) / a;
53 else D = (-2 * cross + a * r2) / (T + pt);
54 //
55 Par(0) = D; // Store D
56 Par(1) = phi0; // Store phi0
57 Par(2) = C; // Store C
58 //Longitudinal parameters
[c5696dd]59 Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));
60 Double_t st = asin(B) / C;
[82db145]61 Double_t ct = p(2) / pt;
62 Double_t z0 = x(2) - ct * st;
63 //
64 Par(3) = z0; // Store z0
65 Par(4) = ct; // Store cot(theta)
66 //
67 return Par;
68}
69// non-static
70TVectorD TrkUtil::XPtoPar(TVector3 x, TVector3 p, Double_t Q)
71{
72 //
73 TVectorD Par(5);
74 Double_t Bz = fBz;
75 Par = XPtoPar(x, p, Q, Bz);
76 //
77 return Par;
78}
79//
80TVector3 TrkUtil::ParToX(TVectorD Par)
81{
82 Double_t D = Par(0);
83 Double_t phi0 = Par(1);
84 Double_t z0 = Par(3);
85 //
86 TVector3 Xval;
[c5696dd]87 Xval(0) = -D * sin(phi0);
88 Xval(1) = D * cos(phi0);
[82db145]89 Xval(2) = z0;
90 //
91 return Xval;
92}
93//
94TVector3 TrkUtil::ParToP(TVectorD Par)
95{
96 if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl;
97 //
98 return ParToP(Par,fBz);
99}
100//
101TVector3 TrkUtil::ParToP(TVectorD Par, Double_t Bz)
102{
103 Double_t C = Par(2);
104 Double_t phi0 = Par(1);
105 Double_t ct = Par(4);
106 //
107 TVector3 Pval;
108 Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C);
[c5696dd]109 Pval(0) = pt * cos(phi0);
110 Pval(1) = pt * sin(phi0);
[82db145]111 Pval(2) = pt * ct;
112 //
113 return Pval;
114}
115//
116Double_t TrkUtil::ParToQ(TVectorD Par)
117{
118 return TMath::Sign(1.0, -Par(2));
119}
120
121//
122// Parameter conversion to ACTS format
123TVectorD TrkUtil::ParToACTS(TVectorD Par)
124{
125 TVectorD pACTS(6); // Return vector
126 //
127 Double_t b = -cSpeed() * fBz / 2.;
128 pACTS(0) = 1000 * Par(0); // D from m to mm
129 pACTS(1) = 1000 * Par(3); // z0 from m to mm
130 pACTS(2) = Par(1); // Phi0 is unchanged
[c5696dd]131 pACTS(3) = atan2(1.0, Par(4)); // Theta in [0, pi] range
132 pACTS(4) = Par(2) / (b * sqrt(1 + Par(4) * Par(4))); // q/p in GeV
[82db145]133 pACTS(5) = 0.0; // Time: currently undefined
134 //
135 return pACTS;
136}
137// Covariance conversion to ACTS format
138TMatrixDSym TrkUtil::CovToACTS(TVectorD Par, TMatrixDSym Cov)
139{
140 TMatrixDSym cACTS(6); cACTS.Zero();
141 Double_t b = -cSpeed() * fBz / 2.;
142 //
143 // Fill derivative matrix
144 TMatrixD A(5, 5); A.Zero();
145 Double_t ct = Par(4); // cot(theta)
146 Double_t C = Par(2); // half curvature
147 A(0, 0) = 1000.; // D-D conversion to mm
148 A(1, 2) = 1.0; // phi0-phi0
[c5696dd]149 A(2, 4) = 1.0 / (sqrt(1.0 + ct * ct) * b); // q/p-C
[82db145]150 A(3, 1) = 1000.; // z0-z0 conversion to mm
151 A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta)
152 A(4, 4) = -C * ct / (b * pow(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta)
153 //
154 TMatrixDSym Cv = Cov;
155 TMatrixD At(5, 5);
156 At.Transpose(A);
157 Cv.Similarity(At);
158 TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv;
159 cACTS(5, 5) = 0.1; // Currently undefined: set to arbitrary value to avoid crashes
160 //
161 return cACTS;
162}
163//
164// Parameter conversion to ILC format
165TVectorD TrkUtil::ParToILC(TVectorD Par)
166{
167 TVectorD pILC(5); // Return vector
168 //
169 pILC(0) = Par(0) * 1.0e3; // d0 in mm
170 pILC(1) = Par(1); // phi0 is unchanged
171 pILC(2) = -2 * Par(2) * 1.0e-3; // w in mm^-1
172 pILC(3) = Par(3) * 1.0e3; // z0 in mm
173 pILC(4) = Par(4); // tan(lambda) = cot(theta)
174 //
175 return pILC;
176}
177// Covariance conversion to ILC format
178TMatrixDSym TrkUtil::CovToILC(TMatrixDSym Cov)
179{
180 TMatrixDSym cILC(5); cILC.Zero();
181 //
182 // Fill derivative matrix
183 TMatrixD A(5, 5); A.Zero();
184 //
185 A(0, 0) = 1.0e3; // D-d0 in mm
186 A(1, 1) = 1.0; // phi0-phi0
187 A(2, 2) = -2.0e-3; // w-C
188 A(3, 3) = 1.0e3; // z0-z0 conversion to mm
189 A(4, 4) = 1.0; // tan(lambda) - cot(theta)
190 //
191 TMatrixDSym Cv = Cov;
192 TMatrixD At(5, 5);
193 At.Transpose(A);
194 Cv.Similarity(At);
195 cILC = Cv;
196 //
197 return cILC;
198}
199//
200// Conversion from meters to mm
201TVectorD TrkUtil::ParToMm(TVectorD Par) // Parameter conversion
202{
203 TVectorD Pmm(5); // Return vector
204 //
205 Pmm(0) = Par(0) * 1.0e3; // d0 in mm
206 Pmm(1) = Par(1); // phi0 is unchanged
207 Pmm(2) = Par(2) * 1.0e-3; // C in mm^-1
208 Pmm(3) = Par(3) * 1.0e3; // z0 in mm
209 Pmm(4) = Par(4); // tan(lambda) = cot(theta) unchanged
210 //
211 return Pmm;
212}
213TMatrixDSym TrkUtil::CovToMm(TMatrixDSym Cov) // Covariance conversion
214{
215 TMatrixDSym Cmm(5); Cmm.Zero();
216 //
217 // Fill derivative matrix
218 TMatrixD A(5, 5); A.Zero();
219 //
220 A(0, 0) = 1.0e3; // D-d0 in mm
221 A(1, 1) = 1.0; // phi0-phi0
222 A(2, 2) = 1.0e-3; // C-C
223 A(3, 3) = 1.0e3; // z0-z0 conversion to mm
224 A(4, 4) = 1.0; // lambda - cot(theta)
225 //
226 TMatrixDSym Cv = Cov;
227 TMatrixD At(5, 5);
228 At.Transpose(A);
229 Cv.Similarity(At);
230 Cmm = Cv;
231 //
232 return Cmm;
233}
234//
235// Setup chamber volume
236void TrkUtil::SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax)
237{
238 fRmin = Rmin; // Lower DCH radius
239 fRmax = Rmax; // Higher DCH radius
240 fZmin = Zmin; // Lower DCH z
241 fZmax = Zmax; // Higher DCH z
242}
243//
244// Get Trakck length inside DCH volume
245Double_t TrkUtil::TrkLen(TVectorD Par)
246{
247 Double_t tLength = 0.0;
248 // Check if geometry is initialized
249 if (fZmin == 0.0 && fZmax == 0.0)
250 {
251 // No geometry set so send a warning and return 0
252 std::cout << "TrkUtil::TrkLen() called without a DCH volume defined" << std::endl;
253 }
254 else
255 {
256 //******************************************************************
257 // Determine the track length inside the chamber ****
258 //******************************************************************
259 //
260 // Track pararameters
261 Double_t D = Par(0); // Transverse impact parameter
262 Double_t phi0 = Par(1); // Transverse direction at minimum approach
263 Double_t C = Par(2); // Half curvature
264 Double_t z0 = Par(3); // Z at minimum approach
265 Double_t ct = Par(4); // cot(theta)
266 //std::cout << "TrkUtil:: parameters: D= " << D << ", phi0= " << phi0
267 // << ", C= " << C << ", z0= " << z0 << ", ct= " << ct << std::endl;
268 //
269 // Track length per unit phase change
[c5696dd]270 Double_t Scale = sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C));
[82db145]271 //
272 // Find intersections with chamber boundaries
273 //
274 Double_t phRin = 0.0; // phase of inner cylinder
275 Double_t phRin2= 0.0; // phase of inner cylinder intersection (2nd branch)
276 Double_t phRhi = 0.0; // phase of outer cylinder intersection
277 Double_t phZmn = 0.0; // phase of left wall intersection
278 Double_t phZmx = 0.0; // phase of right wall intersection
279 // ... with inner cylinder
280 Double_t Rtop = TMath::Abs((1.0 + C*D) / C);
281
282 if (Rtop > fRmin && TMath::Abs(D) < fRmin) // *** don't treat large D tracks for the moment ***
283 {
[c5696dd]284 Double_t ph = 2 * asin(C*sqrt((fRmin*fRmin - D*D) / (1.0 + 2.0*C*D)));
[82db145]285 Double_t z = z0 + ct*ph / (2.0*C);
286
287 //std::cout << "Rin intersection: ph = " << ph<<", z= "<<z << std::endl;
288
289 if (z < fZmax && z > fZmin) phRin = TMath::Abs(ph); // Intersection inside chamber volume
290 //
291 // Include second branch of loopers
[59ba063]292 Double_t Pi = 3.14159265358979323846;
293 Double_t ph2 = 2*Pi - TMath::Abs(ph);
[82db145]294 if (ph < 0)ph2 = -ph2;
295 z = z0 + ct * ph2 / (2.0 * C);
296 if (z < fZmax && z > fZmin) phRin2 = TMath::Abs(ph2); // Intersection inside chamber volume
297 }
298 // ... with outer cylinder
299 if (Rtop > fRmax && TMath::Abs(D) < fRmax) // *** don't treat large D tracks for the moment ***
300 {
[c5696dd]301 Double_t ph = 2 * asin(C*sqrt((fRmax*fRmax - D*D) / (1.0 + 2.0*C*D)));
[82db145]302 Double_t z = z0 + ct*ph / (2.0*C);
303 if (z < fZmax && z > fZmin) phRhi = TMath::Abs(ph); // Intersection inside chamber volume
304 }
305 // ... with left wall
306 Double_t Zdir = (fZmin - z0) / ct;
307 if (Zdir > 0.0)
308 {
309 Double_t ph = 2.0*C*Zdir;
[c5696dd]310 Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C));
[82db145]311 if (Rint < fRmax && Rint > fRmin) phZmn = TMath::Abs(ph); // Intersection inside chamber volume
312 }
313 // ... with right wall
314 Zdir = (fZmax - z0) / ct;
315 if (Zdir > 0.0)
316 {
317 Double_t ph = 2.0*C*Zdir;
[c5696dd]318 Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C));
[82db145]319 if (Rint < fRmax && Rint > fRmin) phZmx = TMath::Abs(ph); // Intersection inside chamber volume
320 }
321 //
322 // Order phases and keep the lowest two non-zero ones
323 //
324 const Int_t Nint = 5;
325 Double_t dPhase = 0.0; // Phase difference between two close intersections
326 Double_t ph_arr[Nint] = { phRin, phRin2, phRhi, phZmn, phZmx };
[59ba063]327 std::sort(ph_arr, ph_arr + Nint);
[82db145]328 Int_t iPos = -1; // First element > 0
329 for (Int_t i = 0; i < Nint; i++)
330 {
[59ba063]331 if (ph_arr[i] <= 0.0) iPos = i;
[82db145]332 }
333
334 if (iPos < Nint - 2)
335 {
[59ba063]336 dPhase = ph_arr[iPos + 2] - ph_arr[iPos + 1];
[82db145]337 tLength = dPhase*Scale;
338 }
339 }
340 return tLength;
341}
342//
343// Return number of ionization clusters
344Bool_t TrkUtil::IonClusters(Double_t &Ncl, Double_t mass, TVectorD Par)
345{
346 //
347 // Units are meters/Tesla/GeV
348 //
349 Ncl = 0.0;
350 Bool_t Signal = kFALSE;
351 Double_t tLen = 0;
352 // Check if geometry is initialized
353 if (fZmin == 0.0 && fZmax == 0.0)
354 {
355 // No geometry set so send a warning and return 0
356 std::cout << "TrkUtil::IonClusters() called without a volume defined" << std::endl;
357 }
358 else tLen = TrkLen(Par);
359
360 //******************************************************************
361 // Now get the number of clusters ****
362 //******************************************************************
363 //
364 Double_t muClu = 0.0; // mean number of clusters
365 Double_t bg = 0.0; // beta*gamma
366 Ncl = 0.0;
367 if (tLen > 0.0)
368 {
369 Signal = kTRUE;
370 //
371 // Find beta*gamma
372 if (fBz == 0.0)
373 {
374 Signal = kFALSE;
375 std::cout << "TrkUtil::IonClusters: Please set Bz!!!" << std::endl;
376 }
377 else
378 {
379 TVector3 p = ParToP(Par);
380 bg = p.Mag() / mass;
381 muClu = Nclusters(bg)*tLen; // Avg. number of clusters
382
383 Ncl = gRandom->PoissonD(muClu); // Actual number of clusters
384 }
385
386 }
387//
388 return Signal;
389}
390//
391//
392Double_t TrkUtil::Nclusters(Double_t begam)
393{
394 Int_t Opt = fGasSel;
395 Double_t Nclu = Nclusters(begam, Opt);
396 //
397 return Nclu;
398}
399//
400Double_t TrkUtil::Nclusters(Double_t begam, Int_t Opt) {
401 //
402 // Opt = 0: He 90 - Isobutane 10
403 // = 1: pure He
404 // = 2: Argon 50 - Ethane 50
405 // = 3: pure Argon
406 //
407 //
[4df491e]408 /*
[82db145]409 std::vector<double> bg{ 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10.,
410 12., 15., 20., 50., 100., 200., 500., 1000. };
411 // He 90 - Isobutane 10
412 std::vector<double> ncl_He_Iso{ 42.94, 23.6,18.97,12.98,12.2,12.13,
413 12.24,12.73,13.03,13.29,13.63,14.08,15.56,16.43,16.8,16.95,16.98 };
414 //
415 // pure He
416 std::vector<double> ncl_He{ 11.79,6.5,5.23,3.59,3.38,3.37,3.4,3.54,3.63,
417 3.7,3.8,3.92,4.33,4.61,4.78,4.87,4.89 };
418 //
419 // Argon 50 - Ethane 50
420 std::vector<double> ncl_Ar_Eth{ 130.04,71.55,57.56,39.44,37.08,36.9,
421 37.25,38.76,39.68,40.49,41.53,42.91,46.8,48.09,48.59,48.85,48.93 };
422 //
423 // pure Argon
424 std::vector<double> ncl_Ar{ 88.69,48.93,39.41,27.09,25.51,25.43,25.69,
425 26.78,27.44,28.02,28.77,29.78,32.67,33.75,34.24,34.57,34.68 };
426 //
427 Int_t nPoints = (Int_t)bg.size();
428 bg.push_back(10000.);
429 std::vector<double> ncl;
430 switch (Opt)
431 {
432 case 0: ncl = ncl_He_Iso; // He-Isobutane
433 break;
434 case 1: ncl = ncl_He; // pure He
435 break;
436 case 2: ncl = ncl_Ar_Eth; // Argon - Ethane
437 break;
438 case 3: ncl = ncl_Ar; // pure Argon
439 break;
440 }
441 ncl.push_back(ncl[nPoints - 1]);
[4df491e]442 */
443 const Int_t Npt = 18;
444 Double_t bg[Npt] = { 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10.,
445 12., 15., 20., 50., 100., 200., 500., 1000., 10000. };
446 //
447 // He 90 - Isobutane 10
448 Double_t ncl_He_Iso[Npt] = { 42.94, 23.6,18.97,12.98,12.2,12.13,
449 12.24,12.73,13.03,13.29,13.63,14.08,15.56,16.43,16.8,16.95,16.98, 16.98 };
450 //
451 // pure He
452 Double_t ncl_He[Npt] = { 11.79,6.5,5.23,3.59,3.38,3.37,3.4,3.54,3.63,
453 3.7,3.8,3.92,4.33,4.61,4.78,4.87,4.89, 4.89 };
454 //
455 // Argon 50 - Ethane 50
456 Double_t ncl_Ar_Eth[Npt] = { 130.04,71.55,57.56,39.44,37.08,36.9,
457 37.25,38.76,39.68,40.49,41.53,42.91,46.8,48.09,48.59,48.85,48.93,48.93 };
458 //
459 // pure Argon
460 Double_t ncl_Ar[Npt] = { 88.69,48.93,39.41,27.09,25.51,25.43,25.69,
461 26.78,27.44,28.02,28.77,29.78,32.67,33.75,34.24,34.57,34.68, 34.68 };
462 //
463 Double_t ncl[Npt];
464 switch (Opt)
465 {
466 case 0: std::copy(ncl_He_Iso, ncl_He_Iso + Npt, ncl); // He-Isobutane
467 break;
468 case 1: std::copy(ncl_He, ncl_He + Npt, ncl); // pure He
469 break;
470 case 2: std::copy(ncl_Ar_Eth, ncl_Ar_Eth + Npt, ncl); // Argon - Ethane
471 break;
472 case 3: std::copy(ncl_Ar, ncl_Ar + Npt, ncl); // pure Argon
473 break;
474 }
475 //
[82db145]476 Int_t ilow = 0;
477 while (begam > bg[ilow])ilow++;
478 ilow--;
479 //std::cout << "ilow= " << ilow << ", low = " << bg[ilow] << ", val = " << begam
480 // << ", high = " << bg[ilow + 1] << std::endl;
481 //
482 Int_t ind[3] = { ilow, ilow + 1, ilow + 2 };
483 TVectorD y(3);
484 for (Int_t i = 0; i < 3; i++)y(i) = ncl[ind[i]];
485 TVectorD x(3);
486 for (Int_t i = 0; i < 3; i++)x(i) = bg[ind[i]];
487 TMatrixD Xval(3, 3);
488 for (Int_t i = 0; i < 3; i++)Xval(i, 0) = 1.0;
489 for (Int_t i = 0; i < 3; i++)Xval(i, 1) = x(i);
490 for (Int_t i = 0; i < 3; i++)Xval(i, 2) = x(i) * x(i);
491 //std::cout << "Xval:" << std::endl; Xval.Print();
492 Xval.Invert();
493 TVectorD coeff = Xval * y;
494 Double_t interp = coeff[0] + coeff[1] * begam + coeff[2] * begam * begam;
495 //std::cout << "val1= (" <<x(0)<<", "<< y(0) << "), val2= ("
496 // <<x(1)<<", "<< y(1) << "), val3= ("
497 // <<x(2)<<", "<< y(2)
498 // << "), result= (" <<begam<<", "<< interp<<")" << std::endl;
499 //
[c5696dd]500 //if (TMath::IsNaN(interp))std::cout << "NaN found: bg= " << begam << ", Opt= " << Opt << std::endl;
[82db145]501 if (begam < bg[0]) interp = 0.0;
502 //std::cout << "bg= " << begam << ", Opt= " << Opt <<", interp = "<<interp<< std::endl;
503 return 100*interp;
504}
505//
506Double_t TrkUtil::funcNcl(Double_t *xp, Double_t *par){
507 Double_t bg = xp[0];
508 return Nclusters(bg);
509}
510//
511void TrkUtil::SetGasMix(Int_t Opt)
512{
513 if (Opt < 0 || Opt > 3)
514 {
515 std::cout << "TrkUtil::SetGasMix Gas option not allowed. No action."
516 << std::endl;
517 }
518 else fGasSel = Opt;
[c5696dd]519}
Note: See TracBrowser for help on using the repository browser.