Changeset 82db145 in git for external/TrackCovariance/TrkUtil.cc
- Timestamp:
- Apr 12, 2021, 6:33:23 PM (4 years ago)
- Branches:
- master
- Children:
- c5696dd
- Parents:
- 3cfe61d
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/TrkUtil.cc
r3cfe61d r82db145 1 1 #include "TrkUtil.h" 2 #include <TMath.h>3 2 #include <iostream> 4 3 … … 7 6 { 8 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 9 13 } 10 14 TrkUtil::TrkUtil() 11 15 { 12 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 13 22 } 14 23 // … … 17 26 { 18 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 19 33 } 20 34 // … … 29 43 Double_t pt = p.Pt(); 30 44 Double_t C = a / (2 * pt); // Half curvature 31 // cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C <<endl;45 //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl; 32 46 Double_t r2 = x.Perp2(); 33 47 Double_t cross = x(0) * p(1) - x(1) * p(0); … … 79 93 TVector3 TrkUtil::ParToP(TVectorD Par) 80 94 { 81 if (fBz == 0.0) 82 std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl; 95 if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl; 83 96 // 84 97 return ParToP(Par,fBz); … … 113 126 Double_t b = -cSpeed() * fBz / 2.; 114 127 pACTS(0) = 1000 * Par(0); // D from m to mm 115 pACTS(1) = 1000 * Par(3); 128 pACTS(1) = 1000 * Par(3); // z0 from m to mm 116 129 pACTS(2) = Par(1); // Phi0 is unchanged 117 130 pACTS(3) = TMath::ATan2(1.0, Par(4)); // Theta in [0, pi] range … … 136 149 A(3, 1) = 1000.; // z0-z0 conversion to mm 137 150 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)151 A(4, 4) = -C * ct / (b * pow(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta) 139 152 // 140 153 TMatrixDSym Cv = Cov; … … 218 231 return Cmm; 219 232 } 233 // 234 // Setup chamber volume 235 void 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 244 Double_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 343 Bool_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 // 391 Double_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 // 399 Double_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 // 470 Double_t TrkUtil::funcNcl(Double_t *xp, Double_t *par){ 471 Double_t bg = xp[0]; 472 return Nclusters(bg); 473 } 474 // 475 void 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 TracChangeset
for help on using the changeset viewer.