Changes in external/TrackCovariance/TrkUtil.cc [f9517a5:4df491e] in git
- File:
-
- 1 edited
-
external/TrackCovariance/TrkUtil.cc (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/TrkUtil.cc
rf9517a5 r4df491e 1 1 #include "TrkUtil.h" 2 #include <TMath.h>3 2 #include <iostream> 3 #include <algorithm> 4 4 5 5 // Constructor … … 7 7 { 8 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 9 14 } 10 15 TrkUtil::TrkUtil() 11 16 { 12 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 13 23 } 14 24 // … … 17 27 { 18 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 19 34 } 20 35 // … … 29 44 Double_t pt = p.Pt(); 30 45 Double_t C = a / (2 * pt); // Half curvature 31 // cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C <<endl;46 //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl; 32 47 Double_t r2 = x.Perp2(); 33 48 Double_t cross = x(0) * p(1) - x(1) * p(0); 34 Double_t T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2);35 Double_t phi0 = TMath::ATan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); // Phi049 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 36 51 Double_t D; // Impact parameter D 37 52 if (pt < 10.0) D = (T - pt) / a; … … 42 57 Par(2) = C; // Store C 43 58 //Longitudinal parameters 44 Double_t B = C * TMath::Sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));45 Double_t st = TMath::ASin(B) / C;59 Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D)); 60 Double_t st = asin(B) / C; 46 61 Double_t ct = p(2) / pt; 47 62 Double_t z0 = x(2) - ct * st; … … 70 85 // 71 86 TVector3 Xval; 72 Xval(0) = -D * TMath::Sin(phi0);73 Xval(1) = D * TMath::Cos(phi0);87 Xval(0) = -D * sin(phi0); 88 Xval(1) = D * cos(phi0); 74 89 Xval(2) = z0; 75 90 // … … 79 94 TVector3 TrkUtil::ParToP(TVectorD Par) 80 95 { 81 if (fBz == 0.0) 82 std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl; 96 if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl; 83 97 // 84 98 return ParToP(Par,fBz); … … 93 107 TVector3 Pval; 94 108 Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C); 95 Pval(0) = pt * TMath::Cos(phi0);96 Pval(1) = pt * TMath::Sin(phi0);109 Pval(0) = pt * cos(phi0); 110 Pval(1) = pt * sin(phi0); 97 111 Pval(2) = pt * ct; 98 112 // … … 113 127 Double_t b = -cSpeed() * fBz / 2.; 114 128 pACTS(0) = 1000 * Par(0); // D from m to mm 115 pACTS(1) = 1000 * Par(3); 129 pACTS(1) = 1000 * Par(3); // z0 from m to mm 116 130 pACTS(2) = Par(1); // Phi0 is unchanged 117 pACTS(3) = TMath::ATan2(1.0, Par(4)); // Theta in [0, pi] range118 pACTS(4) = Par(2) / (b * TMath::Sqrt(1 + Par(4) * Par(4))); // q/p in GeV131 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 119 133 pACTS(5) = 0.0; // Time: currently undefined 120 134 // … … 133 147 A(0, 0) = 1000.; // D-D conversion to mm 134 148 A(1, 2) = 1.0; // phi0-phi0 135 A(2, 4) = 1.0 / ( TMath::Sqrt(1.0 + ct * ct) * b); // q/p-C149 A(2, 4) = 1.0 / (sqrt(1.0 + ct * ct) * b); // q/p-C 136 150 A(3, 1) = 1000.; // z0-z0 conversion to mm 137 151 A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta) 138 A(4, 4) = -C * ct / (b * TMath::Power(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta)152 A(4, 4) = -C * ct / (b * pow(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta) 139 153 // 140 154 TMatrixDSym Cv = Cov; … … 218 232 return Cmm; 219 233 } 234 // 235 // Setup chamber volume 236 void 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 245 Double_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 270 Double_t Scale = sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C)); 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 { 284 Double_t ph = 2 * asin(C*sqrt((fRmin*fRmin - D*D) / (1.0 + 2.0*C*D))); 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 292 Double_t Pi = 3.14159265358979323846; 293 Double_t ph2 = 2*Pi - TMath::Abs(ph); 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 { 301 Double_t ph = 2 * asin(C*sqrt((fRmax*fRmax - D*D) / (1.0 + 2.0*C*D))); 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; 310 Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C)); 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; 318 Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C)); 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 }; 327 std::sort(ph_arr, ph_arr + Nint); 328 Int_t iPos = -1; // First element > 0 329 for (Int_t i = 0; i < Nint; i++) 330 { 331 if (ph_arr[i] <= 0.0) iPos = i; 332 } 333 334 if (iPos < Nint - 2) 335 { 336 dPhase = ph_arr[iPos + 2] - ph_arr[iPos + 1]; 337 tLength = dPhase*Scale; 338 } 339 } 340 return tLength; 341 } 342 // 343 // Return number of ionization clusters 344 Bool_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 // 392 Double_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 // 400 Double_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 // 408 /* 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]); 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 // 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 // 500 //if (TMath::IsNaN(interp))std::cout << "NaN found: bg= " << begam << ", Opt= " << Opt << std::endl; 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 // 506 Double_t TrkUtil::funcNcl(Double_t *xp, Double_t *par){ 507 Double_t bg = xp[0]; 508 return Nclusters(bg); 509 } 510 // 511 void 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; 519 }
Note:
See TracChangeset
for help on using the changeset viewer.
