Changeset 65776c0 in git for external/TrackCovariance/TrkUtil.cc
- Timestamp:
- May 6, 2021, 10:39:31 AM (4 years ago)
- Branches:
- master
- Children:
- 781af69
- Parents:
- 7fbde86
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/TrkUtil.cc
r7fbde86 r65776c0 2 2 #include <iostream> 3 3 #include <algorithm> 4 #include <TSpline.h> 4 5 5 6 // Constructor … … 45 46 Double_t C = a / (2 * pt); // Half curvature 46 47 //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl; 47 Double_t r2 = x .Perp2();48 Double_t r2 = x(0) * x(0) + x(1) * x(1); 48 49 Double_t cross = x(0) * p(1) - x(1) * p(0); 49 50 Double_t T = sqrt(pt * pt - 2 * a * cross + a * a * r2); … … 60 61 Double_t st = asin(B) / C; 61 62 Double_t ct = p(2) / pt; 62 Double_t z0 = x(2) - ct * st; 63 Double_t z0; 64 Double_t dot = x(0) * p(0) + x(1) * p(1); 65 if (dot > 0.0) z0 = x(2) - ct * st; 66 else z0 = x(2) + ct * st; 63 67 // 64 68 Par(3) = z0; // Store z0 … … 96 100 if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl; 97 101 // 98 return ParToP(Par, fBz);102 return ParToP(Par, fBz); 99 103 } 100 104 // … … 233 237 } 234 238 // 235 236 //237 void TrkUtil::SetBfield(Double_t Bz)238 {239 fBz = Bz;240 }241 242 239 // Setup chamber volume 243 240 void TrkUtil::SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax) … … 274 271 // << ", C= " << C << ", z0= " << z0 << ", ct= " << ct << std::endl; 275 272 // 276 // Track length per unit phase change 277 Double_t Scale = sqrt(1.0 + ct *ct) / (2.0*TMath::Abs(C));273 // Track length per unit phase change 274 Double_t Scale = sqrt(1.0 + ct * ct) / (2.0 * TMath::Abs(C)); 278 275 // 279 276 // Find intersections with chamber boundaries 280 277 // 281 Double_t phRin = 0.0; // phase of inner cylinder 282 Double_t phRin2 = 0.0; // phase of inner cylinder intersection (2nd branch)278 Double_t phRin = 0.0; // phase of inner cylinder 279 Double_t phRin2 = 0.0; // phase of inner cylinder intersection (2nd branch) 283 280 Double_t phRhi = 0.0; // phase of outer cylinder intersection 284 281 Double_t phZmn = 0.0; // phase of left wall intersection 285 282 Double_t phZmx = 0.0; // phase of right wall intersection 286 283 // ... with inner cylinder 287 Double_t Rtop = TMath::Abs((1.0 + C *D) / C);284 Double_t Rtop = TMath::Abs((1.0 + C * D) / C); 288 285 289 286 if (Rtop > fRmin && TMath::Abs(D) < fRmin) // *** don't treat large D tracks for the moment *** 290 287 { 291 Double_t ph = 2 * asin(C *sqrt((fRmin*fRmin - D*D) / (1.0 + 2.0*C*D)));292 Double_t z = z0 + ct *ph / (2.0*C);288 Double_t ph = 2 * asin(C * sqrt((fRmin * fRmin - D * D) / (1.0 + 2.0 * C * D))); 289 Double_t z = z0 + ct * ph / (2.0 * C); 293 290 294 291 //std::cout << "Rin intersection: ph = " << ph<<", z= "<<z << std::endl; 295 292 296 if (z < fZmax && z > fZmin) phRin = TMath::Abs(ph); // Intersection inside chamber volume 293 if (z < fZmax && z > fZmin) phRin = TMath::Abs(ph); // Intersection inside chamber volume 297 294 // 298 295 // Include second branch of loopers 299 296 Double_t Pi = 3.14159265358979323846; 300 Double_t ph2 = 2 *Pi - TMath::Abs(ph);297 Double_t ph2 = 2 * Pi - TMath::Abs(ph); 301 298 if (ph < 0)ph2 = -ph2; 302 299 z = z0 + ct * ph2 / (2.0 * C); … … 306 303 if (Rtop > fRmax && TMath::Abs(D) < fRmax) // *** don't treat large D tracks for the moment *** 307 304 { 308 Double_t ph = 2 * asin(C *sqrt((fRmax*fRmax - D*D) / (1.0 + 2.0*C*D)));309 Double_t z = z0 + ct *ph / (2.0*C);310 if (z < fZmax && z > fZmin) phRhi = TMath::Abs(ph); // Intersection inside chamber volume 305 Double_t ph = 2 * asin(C * sqrt((fRmax * fRmax - D * D) / (1.0 + 2.0 * C * D))); 306 Double_t z = z0 + ct * ph / (2.0 * C); 307 if (z < fZmax && z > fZmin) phRhi = TMath::Abs(ph); // Intersection inside chamber volume 311 308 } 312 309 // ... with left wall … … 314 311 if (Zdir > 0.0) 315 312 { 316 Double_t ph = 2.0 *C*Zdir;317 Double_t Rint = sqrt(D *D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C));318 if (Rint < fRmax && Rint > fRmin) phZmn = TMath::Abs(ph); // Intersection inside chamber volume 313 Double_t ph = 2.0 * C * Zdir; 314 Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C)); 315 if (Rint < fRmax && Rint > fRmin) phZmn = TMath::Abs(ph); // Intersection inside chamber volume 319 316 } 320 317 // ... with right wall … … 322 319 if (Zdir > 0.0) 323 320 { 324 Double_t ph = 2.0 *C*Zdir;325 Double_t Rint = sqrt(D *D + (1.0 + 2.0*C*D)*pow(sin(ph / 2), 2) / (C*C));326 if (Rint < fRmax && Rint > fRmin) phZmx = TMath::Abs(ph); // Intersection inside chamber volume 321 Double_t ph = 2.0 * C * Zdir; 322 Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C)); 323 if (Rint < fRmax && Rint > fRmin) phZmx = TMath::Abs(ph); // Intersection inside chamber volume 327 324 } 328 325 // … … 342 339 { 343 340 dPhase = ph_arr[iPos + 2] - ph_arr[iPos + 1]; 344 tLength = dPhase *Scale;341 tLength = dPhase * Scale; 345 342 } 346 343 } … … 349 346 // 350 347 // Return number of ionization clusters 351 Bool_t TrkUtil::IonClusters(Double_t &Ncl, Double_t mass, TVectorD Par)348 Bool_t TrkUtil::IonClusters(Double_t& Ncl, Double_t mass, TVectorD Par) 352 349 { 353 350 // … … 386 383 TVector3 p = ParToP(Par); 387 384 bg = p.Mag() / mass; 388 muClu = Nclusters(bg) *tLen; // Avg. number of clusters385 muClu = Nclusters(bg) * tLen; // Avg. number of clusters 389 386 390 387 Ncl = gRandom->PoissonD(muClu); // Actual number of clusters … … 392 389 393 390 } 394 //391 // 395 392 return Signal; 396 393 } … … 413 410 // 414 411 // 415 /*416 std::vector<double> bg{ 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10.,417 12., 15., 20., 50., 100., 200., 500., 1000. };418 // He 90 - Isobutane 10419 std::vector<double> ncl_He_Iso{ 42.94, 23.6,18.97,12.98,12.2,12.13,420 12.24,12.73,13.03,13.29,13.63,14.08,15.56,16.43,16.8,16.95,16.98 };421 //422 // pure He423 std::vector<double> ncl_He{ 11.79,6.5,5.23,3.59,3.38,3.37,3.4,3.54,3.63,424 3.7,3.8,3.92,4.33,4.61,4.78,4.87,4.89 };425 //426 // Argon 50 - Ethane 50427 std::vector<double> ncl_Ar_Eth{ 130.04,71.55,57.56,39.44,37.08,36.9,428 37.25,38.76,39.68,40.49,41.53,42.91,46.8,48.09,48.59,48.85,48.93 };429 //430 // pure Argon431 std::vector<double> ncl_Ar{ 88.69,48.93,39.41,27.09,25.51,25.43,25.69,432 26.78,27.44,28.02,28.77,29.78,32.67,33.75,34.24,34.57,34.68 };433 //434 Int_t nPoints = (Int_t)bg.size();435 bg.push_back(10000.);436 std::vector<double> ncl;437 switch (Opt)438 {439 case 0: ncl = ncl_He_Iso; // He-Isobutane440 break;441 case 1: ncl = ncl_He; // pure He442 break;443 case 2: ncl = ncl_Ar_Eth; // Argon - Ethane444 break;445 case 3: ncl = ncl_Ar; // pure Argon446 break;447 }448 ncl.push_back(ncl[nPoints - 1]);449 */450 412 const Int_t Npt = 18; 451 413 Double_t bg[Npt] = { 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10., … … 469 431 // 470 432 Double_t ncl[Npt]; 471 472 473 433 switch (Opt) 434 { 435 case 0: std::copy(ncl_He_Iso, ncl_He_Iso + Npt, ncl); // He-Isobutane 474 436 break; 475 437 case 1: std::copy(ncl_He, ncl_He + Npt, ncl); // pure He 476 438 break; 477 439 case 2: std::copy(ncl_Ar_Eth, ncl_Ar_Eth + Npt, ncl); // Argon - Ethane 478 440 break; 479 441 case 3: std::copy(ncl_Ar, ncl_Ar + Npt, ncl); // pure Argon 480 442 break; 481 } 482 // 483 Int_t ilow = 0; 484 while (begam > bg[ilow])ilow++; 485 ilow--; 486 //std::cout << "ilow= " << ilow << ", low = " << bg[ilow] << ", val = " << begam 487 // << ", high = " << bg[ilow + 1] << std::endl; 488 // 489 Int_t ind[3] = { ilow, ilow + 1, ilow + 2 }; 490 TVectorD y(3); 491 for (Int_t i = 0; i < 3; i++)y(i) = ncl[ind[i]]; 492 TVectorD x(3); 493 for (Int_t i = 0; i < 3; i++)x(i) = bg[ind[i]]; 494 TMatrixD Xval(3, 3); 495 for (Int_t i = 0; i < 3; i++)Xval(i, 0) = 1.0; 496 for (Int_t i = 0; i < 3; i++)Xval(i, 1) = x(i); 497 for (Int_t i = 0; i < 3; i++)Xval(i, 2) = x(i) * x(i); 498 //std::cout << "Xval:" << std::endl; Xval.Print(); 499 Xval.Invert(); 500 TVectorD coeff = Xval * y; 501 Double_t interp = coeff[0] + coeff[1] * begam + coeff[2] * begam * begam; 502 //std::cout << "val1= (" <<x(0)<<", "<< y(0) << "), val2= (" 503 // <<x(1)<<", "<< y(1) << "), val3= (" 504 // <<x(2)<<", "<< y(2) 505 // << "), result= (" <<begam<<", "<< interp<<")" << std::endl; 506 // 507 //if (TMath::IsNaN(interp))std::cout << "NaN found: bg= " << begam << ", Opt= " << Opt << std::endl; 508 if (begam < bg[0]) interp = 0.0; 509 //std::cout << "bg= " << begam << ", Opt= " << Opt <<", interp = "<<interp<< std::endl; 510 return 100*interp; 511 } 512 // 513 Double_t TrkUtil::funcNcl(Double_t *xp, Double_t *par){ 443 } 444 // 445 Double_t interp = 0.0; 446 TSpline3* sp3 = new TSpline3("sp3", bg, ncl, Npt); 447 if (begam > bg[0] && begam < bg[Npt - 1]) interp = sp3->Eval(begam); 448 return 100 * interp; 449 } 450 // 451 Double_t TrkUtil::funcNcl(Double_t* xp, Double_t* par) { 514 452 Double_t bg = xp[0]; 515 453 return Nclusters(bg);
Note:
See TracChangeset
for help on using the changeset viewer.