Changes in external/TrackCovariance/TrkUtil.cc [4df491e:f9517a5] in git
- File:
-
- 1 edited
-
external/TrackCovariance/TrkUtil.cc (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/TrkUtil.cc
r4df491e rf9517a5 1 1 #include "TrkUtil.h" 2 #include <TMath.h> 2 3 #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 radius11 fRmax = 0.0; // Higher DCH radius12 fZmin = 0.0; // Lower DCH z13 fZmax = 0.0; // Higher DCH z14 9 } 15 10 TrkUtil::TrkUtil() 16 11 { 17 12 fBz = 0.0; 18 fGasSel = 0; // Default is He-Isobuthane (90-10)19 fRmin = 0.0; // Lower DCH radius20 fRmax = 0.0; // Higher DCH radius21 fZmin = 0.0; // Lower DCH z22 fZmax = 0.0; // Higher DCH z23 13 } 24 14 // … … 27 17 { 28 18 fBz = 0.0; 29 fGasSel = 0; // Default is He-Isobuthane (90-10)30 fRmin = 0.0; // Lower DCH radius31 fRmax = 0.0; // Higher DCH radius32 fZmin = 0.0; // Lower DCH z33 fZmax = 0.0; // Higher DCH z34 19 } 35 20 // … … 44 29 Double_t pt = p.Pt(); 45 30 Double_t C = a / (2 * pt); // Half curvature 46 // std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl;31 //cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << endl; 47 32 Double_t r2 = x.Perp2(); 48 33 Double_t cross = x(0) * p(1) - x(1) * p(0); 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); // Phi034 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); // Phi0 51 36 Double_t D; // Impact parameter D 52 37 if (pt < 10.0) D = (T - pt) / a; … … 57 42 Par(2) = C; // Store C 58 43 //Longitudinal parameters 59 Double_t B = C * sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D));60 Double_t st = asin(B) / C;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; 61 46 Double_t ct = p(2) / pt; 62 47 Double_t z0 = x(2) - ct * st; … … 85 70 // 86 71 TVector3 Xval; 87 Xval(0) = -D * sin(phi0);88 Xval(1) = D * cos(phi0);72 Xval(0) = -D * TMath::Sin(phi0); 73 Xval(1) = D * TMath::Cos(phi0); 89 74 Xval(2) = z0; 90 75 // … … 94 79 TVector3 TrkUtil::ParToP(TVectorD Par) 95 80 { 96 if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl; 81 if (fBz == 0.0) 82 std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl; 97 83 // 98 84 return ParToP(Par,fBz); … … 107 93 TVector3 Pval; 108 94 Double_t pt = Bz * cSpeed() / TMath::Abs(2 * C); 109 Pval(0) = pt * cos(phi0);110 Pval(1) = pt * sin(phi0);95 Pval(0) = pt * TMath::Cos(phi0); 96 Pval(1) = pt * TMath::Sin(phi0); 111 97 Pval(2) = pt * ct; 112 98 // … … 127 113 Double_t b = -cSpeed() * fBz / 2.; 128 114 pACTS(0) = 1000 * Par(0); // D from m to mm 129 pACTS(1) = 1000 * Par(3); // z0 from m to mm115 pACTS(1) = 1000 * Par(3); // z0 from m to mm 130 116 pACTS(2) = Par(1); // Phi0 is unchanged 131 pACTS(3) = atan2(1.0, Par(4)); // Theta in [0, pi] range132 pACTS(4) = Par(2) / (b * sqrt(1 + Par(4) * Par(4))); // q/p in GeV117 pACTS(3) = TMath::ATan2(1.0, Par(4)); // Theta in [0, pi] range 118 pACTS(4) = Par(2) / (b * TMath::Sqrt(1 + Par(4) * Par(4))); // q/p in GeV 133 119 pACTS(5) = 0.0; // Time: currently undefined 134 120 // … … 147 133 A(0, 0) = 1000.; // D-D conversion to mm 148 134 A(1, 2) = 1.0; // phi0-phi0 149 A(2, 4) = 1.0 / ( sqrt(1.0 + ct * ct) * b); // q/p-C135 A(2, 4) = 1.0 / (TMath::Sqrt(1.0 + ct * ct) * b); // q/p-C 150 136 A(3, 1) = 1000.; // z0-z0 conversion to mm 151 137 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)138 A(4, 4) = -C * ct / (b * TMath::Power(1.0 + ct * ct, 3.0 / 2.0)); // q/p-cot(theta) 153 139 // 154 140 TMatrixDSym Cv = Cov; … … 232 218 return Cmm; 233 219 } 234 //235 // Setup chamber volume236 void TrkUtil::SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax)237 {238 fRmin = Rmin; // Lower DCH radius239 fRmax = Rmax; // Higher DCH radius240 fZmin = Zmin; // Lower DCH z241 fZmax = Zmax; // Higher DCH z242 }243 //244 // Get Trakck length inside DCH volume245 Double_t TrkUtil::TrkLen(TVectorD Par)246 {247 Double_t tLength = 0.0;248 // Check if geometry is initialized249 if (fZmin == 0.0 && fZmax == 0.0)250 {251 // No geometry set so send a warning and return 0252 std::cout << "TrkUtil::TrkLen() called without a DCH volume defined" << std::endl;253 }254 else255 {256 //******************************************************************257 // Determine the track length inside the chamber ****258 //******************************************************************259 //260 // Track pararameters261 Double_t D = Par(0); // Transverse impact parameter262 Double_t phi0 = Par(1); // Transverse direction at minimum approach263 Double_t C = Par(2); // Half curvature264 Double_t z0 = Par(3); // Z at minimum approach265 Double_t ct = Par(4); // cot(theta)266 //std::cout << "TrkUtil:: parameters: D= " << D << ", phi0= " << phi0267 // << ", C= " << C << ", z0= " << z0 << ", ct= " << ct << std::endl;268 //269 // Track length per unit phase change270 Double_t Scale = sqrt(1.0 + ct*ct) / (2.0*TMath::Abs(C));271 //272 // Find intersections with chamber boundaries273 //274 Double_t phRin = 0.0; // phase of inner cylinder275 Double_t phRin2= 0.0; // phase of inner cylinder intersection (2nd branch)276 Double_t phRhi = 0.0; // phase of outer cylinder intersection277 Double_t phZmn = 0.0; // phase of left wall intersection278 Double_t phZmx = 0.0; // phase of right wall intersection279 // ... with inner cylinder280 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 volume290 //291 // Include second branch of loopers292 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 volume297 }298 // ... with outer cylinder299 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 volume304 }305 // ... with left wall306 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 volume312 }313 // ... with right wall314 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 volume320 }321 //322 // Order phases and keep the lowest two non-zero ones323 //324 const Int_t Nint = 5;325 Double_t dPhase = 0.0; // Phase difference between two close intersections326 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 > 0329 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 clusters344 Bool_t TrkUtil::IonClusters(Double_t &Ncl, Double_t mass, TVectorD Par)345 {346 //347 // Units are meters/Tesla/GeV348 //349 Ncl = 0.0;350 Bool_t Signal = kFALSE;351 Double_t tLen = 0;352 // Check if geometry is initialized353 if (fZmin == 0.0 && fZmax == 0.0)354 {355 // No geometry set so send a warning and return 0356 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 clusters365 Double_t bg = 0.0; // beta*gamma366 Ncl = 0.0;367 if (tLen > 0.0)368 {369 Signal = kTRUE;370 //371 // Find beta*gamma372 if (fBz == 0.0)373 {374 Signal = kFALSE;375 std::cout << "TrkUtil::IonClusters: Please set Bz!!!" << std::endl;376 }377 else378 {379 TVector3 p = ParToP(Par);380 bg = p.Mag() / mass;381 muClu = Nclusters(bg)*tLen; // Avg. number of clusters382 383 Ncl = gRandom->PoissonD(muClu); // Actual number of clusters384 }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 10403 // = 1: pure He404 // = 2: Argon 50 - Ethane 50405 // = 3: pure Argon406 //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 10412 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 He416 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 50420 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 Argon424 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-Isobutane433 break;434 case 1: ncl = ncl_He; // pure He435 break;436 case 2: ncl = ncl_Ar_Eth; // Argon - Ethane437 break;438 case 3: ncl = ncl_Ar; // pure Argon439 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 10448 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 He452 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 50456 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 Argon460 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-Isobutane467 break;468 case 1: std::copy(ncl_He, ncl_He + Npt, ncl); // pure He469 break;470 case 2: std::copy(ncl_Ar_Eth, ncl_Ar_Eth + Npt, ncl); // Argon - Ethane471 break;472 case 3: std::copy(ncl_Ar, ncl_Ar + Npt, ncl); // pure Argon473 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 = " << begam480 // << ", 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.
