Fork me on GitHub

source: git/external/TrackCovariance/TrkUtil.cc@ 82db145

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

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