- Timestamp:
- Apr 13, 2021, 12:40:06 PM (4 years ago)
- Branches:
- master
- Children:
- 59ba063
- Parents:
- 82db145
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/TrkUtil.cc
r82db145 rc5696dd 1 1 #include "TrkUtil.h" 2 2 #include <iostream> 3 #include <algorithm> 3 4 4 5 // Constructor … … 46 47 Double_t r2 = x.Perp2(); 47 48 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); // 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 50 51 Double_t D; // Impact parameter D 51 52 if (pt < 10.0) D = (T - pt) / a; … … 56 57 Par(2) = C; // Store C 57 58 //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;59 Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D)); 60 Double_t st = asin(B) / C; 60 61 Double_t ct = p(2) / pt; 61 62 Double_t z0 = x(2) - ct * st; … … 84 85 // 85 86 TVector3 Xval; 86 Xval(0) = -D * TMath::Sin(phi0);87 Xval(1) = D * TMath::Cos(phi0);87 Xval(0) = -D * sin(phi0); 88 Xval(1) = D * cos(phi0); 88 89 Xval(2) = z0; 89 90 // … … 106 107 TVector3 Pval; 107 108 Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C); 108 Pval(0) = pt * TMath::Cos(phi0);109 Pval(1) = pt * TMath::Sin(phi0);109 Pval(0) = pt * cos(phi0); 110 Pval(1) = pt * sin(phi0); 110 111 Pval(2) = pt * ct; 111 112 // … … 128 129 pACTS(1) = 1000 * Par(3); // z0 from m to mm 129 130 pACTS(2) = Par(1); // Phi0 is unchanged 130 pACTS(3) = TMath::ATan2(1.0, Par(4)); // Theta in [0, pi] range131 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 132 133 pACTS(5) = 0.0; // Time: currently undefined 133 134 // … … 146 147 A(0, 0) = 1000.; // D-D conversion to mm 147 148 A(1, 2) = 1.0; // phi0-phi0 148 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 149 150 A(3, 1) = 1000.; // z0-z0 conversion to mm 150 151 A(4, 3) = -1.0 / (1.0 + ct * ct); // theta - cot(theta) … … 267 268 // 268 269 // Track length per unit phase change 269 Double_t Scale = TMath::Sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C));270 Double_t Scale = sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C)); 270 271 // 271 272 // Find intersections with chamber boundaries … … 281 282 if (Rtop > fRmin && TMath::Abs(D) < fRmin) // *** don't treat large D tracks for the moment *** 282 283 { 283 Double_t ph = 2 * TMath::ASin(C*TMath::Sqrt((fRmin*fRmin - D*D) / (1.0 + 2.0*C*D)));284 Double_t ph = 2 * asin(C*sqrt((fRmin*fRmin - D*D) / (1.0 + 2.0*C*D))); 284 285 Double_t z = z0 + ct*ph / (2.0*C); 285 286 … … 289 290 // 290 291 // Include second branch of loopers 291 Double_t ph2 = TMath::TwoPi() - TMath::Abs(ph);292 Double_t ph2 = 2*TMath::Pi() - TMath::Abs(ph); 292 293 if (ph < 0)ph2 = -ph2; 293 294 z = z0 + ct * ph2 / (2.0 * C); … … 297 298 if (Rtop > fRmax && TMath::Abs(D) < fRmax) // *** don't treat large D tracks for the moment *** 298 299 { 299 Double_t ph = 2 * TMath::ASin(C*TMath::Sqrt((fRmax*fRmax - D*D) / (1.0 + 2.0*C*D)));300 Double_t ph = 2 * asin(C*sqrt((fRmax*fRmax - D*D) / (1.0 + 2.0*C*D))); 300 301 Double_t z = z0 + ct*ph / (2.0*C); 301 302 if (z < fZmax && z > fZmin) phRhi = TMath::Abs(ph); // Intersection inside chamber volume … … 306 307 { 307 308 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 Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C)); 309 310 if (Rint < fRmax && Rint > fRmin) phZmn = TMath::Abs(ph); // Intersection inside chamber volume 310 311 } … … 314 315 { 315 316 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 Double_t Rint = sqrt(D*D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C)); 317 318 if (Rint < fRmax && Rint > fRmin) phZmx = TMath::Abs(ph); // Intersection inside chamber volume 318 319 } … … 324 325 Double_t ph_arr[Nint] = { phRin, phRin2, phRhi, phZmn, phZmx }; 325 326 Int_t srtind[Nint]; 326 TMath::Sort(Nint, ph_arr, srtind, kFALSE); 327 //TMath::Sort(Nint, ph_arr, srtind, kFALSE); 328 for (Int_t i = 0; i < Nint; i++) { srtind[i] = i; } 329 std::sort(srtind, srtind + Nint, CompareAsc<const Double_t*>(ph_arr)); 327 330 Int_t iPos = -1; // First element > 0 328 331 for (Int_t i = 0; i < Nint; i++) … … 462 465 // << "), result= (" <<begam<<", "<< interp<<")" << std::endl; 463 466 // 464 if (TMath::IsNaN(interp))std::cout << "NaN found: bg= " << begam << ", Opt= " << Opt << std::endl;467 //if (TMath::IsNaN(interp))std::cout << "NaN found: bg= " << begam << ", Opt= " << Opt << std::endl; 465 468 if (begam < bg[0]) interp = 0.0; 466 469 //std::cout << "bg= " << begam << ", Opt= " << Opt <<", interp = "<<interp<< std::endl;
Note:
See TracChangeset
for help on using the changeset viewer.