- Timestamp:
- Apr 30, 2021, 4:30:36 PM (4 years ago)
- Branches:
- master
- Children:
- a47edcc
- Parents:
- 952bbbc
- Location:
- external/TrackCovariance
- Files:
-
- 5 edited
-
ObsTrk.cc (modified) (2 diffs)
-
TrkUtil.cc (modified) (13 diffs)
-
TrkUtil.h (modified) (2 diffs)
-
VertexFit.cc (modified) (21 diffs)
-
VertexFit.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/ObsTrk.cc
r952bbbc r53f4746 14 14 ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC) 15 15 { 16 SetB field(B);16 SetB(B); 17 17 fGC = GC; 18 18 fGenX = x; … … 56 56 ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC) 57 57 { 58 SetB field(B);58 SetB(B); 59 59 fGC = GC; 60 60 fGenX.SetXYZ(x[0],x[1],x[2]); -
external/TrackCovariance/TrkUtil.cc
r952bbbc r53f4746 2 2 #include <iostream> 3 3 #include <algorithm> 4 #include <TSpline.h> 4 5 5 6 // Constructor … … 96 97 if (fBz == 0.0)std::cout << "TrkUtil::ParToP: Warning Bz not set" << std::endl; 97 98 // 98 return ParToP(Par, fBz);99 return ParToP(Par, fBz); 99 100 } 100 101 // … … 233 234 } 234 235 // 235 236 //237 void TrkUtil::SetBfield(Double_t Bz)238 {239 fBz = Bz;240 }241 242 236 // Setup chamber volume 243 237 void TrkUtil::SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax) … … 274 268 // << ", C= " << C << ", z0= " << z0 << ", ct= " << ct << std::endl; 275 269 // 276 // Track length per unit phase change 277 Double_t Scale = sqrt(1.0 + ct *ct) / (2.0*TMath::Abs(C));270 // Track length per unit phase change 271 Double_t Scale = sqrt(1.0 + ct * ct) / (2.0 * TMath::Abs(C)); 278 272 // 279 273 // Find intersections with chamber boundaries 280 274 // 281 Double_t phRin = 0.0; // phase of inner cylinder 282 Double_t phRin2 = 0.0; // phase of inner cylinder intersection (2nd branch)275 Double_t phRin = 0.0; // phase of inner cylinder 276 Double_t phRin2 = 0.0; // phase of inner cylinder intersection (2nd branch) 283 277 Double_t phRhi = 0.0; // phase of outer cylinder intersection 284 278 Double_t phZmn = 0.0; // phase of left wall intersection 285 279 Double_t phZmx = 0.0; // phase of right wall intersection 286 280 // ... with inner cylinder 287 Double_t Rtop = TMath::Abs((1.0 + C *D) / C);281 Double_t Rtop = TMath::Abs((1.0 + C * D) / C); 288 282 289 283 if (Rtop > fRmin && TMath::Abs(D) < fRmin) // *** don't treat large D tracks for the moment *** 290 284 { 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);285 Double_t ph = 2 * asin(C * sqrt((fRmin * fRmin - D * D) / (1.0 + 2.0 * C * D))); 286 Double_t z = z0 + ct * ph / (2.0 * C); 293 287 294 288 //std::cout << "Rin intersection: ph = " << ph<<", z= "<<z << std::endl; 295 289 296 if (z < fZmax && z > fZmin) phRin = TMath::Abs(ph); // Intersection inside chamber volume 290 if (z < fZmax && z > fZmin) phRin = TMath::Abs(ph); // Intersection inside chamber volume 297 291 // 298 292 // Include second branch of loopers 299 293 Double_t Pi = 3.14159265358979323846; 300 Double_t ph2 = 2 *Pi - TMath::Abs(ph);294 Double_t ph2 = 2 * Pi - TMath::Abs(ph); 301 295 if (ph < 0)ph2 = -ph2; 302 296 z = z0 + ct * ph2 / (2.0 * C); … … 306 300 if (Rtop > fRmax && TMath::Abs(D) < fRmax) // *** don't treat large D tracks for the moment *** 307 301 { 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 302 Double_t ph = 2 * asin(C * sqrt((fRmax * fRmax - D * D) / (1.0 + 2.0 * C * D))); 303 Double_t z = z0 + ct * ph / (2.0 * C); 304 if (z < fZmax && z > fZmin) phRhi = TMath::Abs(ph); // Intersection inside chamber volume 311 305 } 312 306 // ... with left wall … … 314 308 if (Zdir > 0.0) 315 309 { 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 310 Double_t ph = 2.0 * C * Zdir; 311 Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C)); 312 if (Rint < fRmax && Rint > fRmin) phZmn = TMath::Abs(ph); // Intersection inside chamber volume 319 313 } 320 314 // ... with right wall … … 322 316 if (Zdir > 0.0) 323 317 { 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 318 Double_t ph = 2.0 * C * Zdir; 319 Double_t Rint = sqrt(D * D + (1.0 + 2.0 * C * D) * pow(sin(ph / 2), 2) / (C * C)); 320 if (Rint < fRmax && Rint > fRmin) phZmx = TMath::Abs(ph); // Intersection inside chamber volume 327 321 } 328 322 // … … 342 336 { 343 337 dPhase = ph_arr[iPos + 2] - ph_arr[iPos + 1]; 344 tLength = dPhase *Scale;338 tLength = dPhase * Scale; 345 339 } 346 340 } … … 349 343 // 350 344 // Return number of ionization clusters 351 Bool_t TrkUtil::IonClusters(Double_t &Ncl, Double_t mass, TVectorD Par)345 Bool_t TrkUtil::IonClusters(Double_t& Ncl, Double_t mass, TVectorD Par) 352 346 { 353 347 // … … 386 380 TVector3 p = ParToP(Par); 387 381 bg = p.Mag() / mass; 388 muClu = Nclusters(bg) *tLen; // Avg. number of clusters382 muClu = Nclusters(bg) * tLen; // Avg. number of clusters 389 383 390 384 Ncl = gRandom->PoissonD(muClu); // Actual number of clusters … … 392 386 393 387 } 394 //388 // 395 389 return Signal; 396 390 } … … 413 407 // 414 408 // 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 409 const Int_t Npt = 18; 451 410 Double_t bg[Npt] = { 0.5, 0.8, 1., 2., 3., 4., 5., 8., 10., … … 469 428 // 470 429 Double_t ncl[Npt]; 471 472 473 430 switch (Opt) 431 { 432 case 0: std::copy(ncl_He_Iso, ncl_He_Iso + Npt, ncl); // He-Isobutane 474 433 break; 475 434 case 1: std::copy(ncl_He, ncl_He + Npt, ncl); // pure He 476 435 break; 477 436 case 2: std::copy(ncl_Ar_Eth, ncl_Ar_Eth + Npt, ncl); // Argon - Ethane 478 437 break; 479 438 case 3: std::copy(ncl_Ar, ncl_Ar + Npt, ncl); // pure Argon 480 439 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){ 440 } 441 // 442 Double_t interp = 0.0; 443 TSpline3* sp3 = new TSpline3("sp3", bg, ncl, Npt); 444 if (begam > bg[0] && begam < bg[Npt - 1]) interp = sp3->Eval(begam); 445 if(begam < bg[0]) interp = bg[0]; 446 if(begam > bg[Npt-1]) interp = bg[Npt-1]; 447 return 100 * interp; 448 } 449 // 450 Double_t TrkUtil::funcNcl(Double_t* xp, Double_t* par) { 514 451 Double_t bg = xp[0]; 515 452 return Nclusters(bg); -
external/TrackCovariance/TrkUtil.h
r952bbbc r53f4746 25 25 // Service routines 26 26 // 27 void SetB(Double_t Bz){ fBz = Bz; } 27 28 TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q); 28 29 TVector3 ParToP(TVectorD Par); … … 70 71 // Cluster counting in gas 71 72 // 72 void SetBfield(Double_t Bz) ;73 void SetBfield(Double_t Bz){ fBz = Bz; } 73 74 // Define gas volume (units = meters) 74 75 void SetDchBoundaries(Double_t Rmin, Double_t Rmax, Double_t Zmin, Double_t Zmax); -
external/TrackCovariance/VertexFit.cc
r952bbbc r53f4746 18 18 fxCst.ResizeTo(3); 19 19 fCovCst.ResizeTo(3, 3); 20 fCovCstInv.ResizeTo(3, 3); 20 21 fXv.ResizeTo(3); 21 22 fcovXv.ResizeTo(3, 3); … … 31 32 fxCst.ResizeTo(3); 32 33 fCovCst.ResizeTo(3, 3); 34 fCovCstInv.ResizeTo(3, 3); 33 35 fXv.ResizeTo(3); 34 36 fcovXv.ResizeTo(3, 3); … … 38 40 TVectorD pr = *trkPar[i]; 39 41 fPar.push_back(new TVectorD(pr)); 42 fParNew.push_back(new TVectorD(pr)); 40 43 TMatrixDSym cv = *trkCov[i]; 41 44 fCov.push_back(new TMatrixDSym(cv)); 45 fCovNew.push_back(new TMatrixDSym(cv)); 42 46 } 43 47 fChi2List.ResizeTo(fNtr); … … 54 58 fxCst.ResizeTo(3); 55 59 fCovCst.ResizeTo(3, 3); 60 fCovCstInv.ResizeTo(3, 3); 56 61 fXv.ResizeTo(3); 57 62 fcovXv.ResizeTo(3, 3); … … 61 66 { 62 67 fPar.push_back(new TVectorD(track[i]->GetObsPar())); 68 fParNew.push_back(new TVectorD(track[i]->GetObsPar())); 63 69 fCov.push_back(new TMatrixDSym(track[i]->GetCov())); 70 fCovNew.push_back(new TMatrixDSym(track[i]->GetCov())); 64 71 } 65 72 } … … 73 80 { 74 81 if (fx0i[i]) { fx0i[i]->Clear(); delete fx0i[i]; } 75 if (fai[i]) { fai[i]->Clear(); delete fai[i]; } 76 if (fDi[i]) { fDi[i]->Clear(); delete fDi[i]; } 77 if (fWi[i]) { fWi[i]->Clear(); delete fWi[i]; } 82 if (fai[i]) { fai[i]->Clear(); delete fai[i]; } 83 if (fdi[i]) { fdi[i]->Clear(); delete fdi[i]; } 84 if (fAti[i]) { fAti[i]->Clear(); delete fAti[i]; } 85 if (fDi[i]) { fDi[i]->Clear(); delete fDi[i]; } 86 if (fWi[i]) { fWi[i]->Clear(); delete fWi[i]; } 78 87 if (fWinvi[i]){ fWinvi[i]->Clear(); delete fWinvi[i]; } 79 88 } … … 81 90 fx0i.clear(); 82 91 fai.clear(); 92 fdi.clear(); 93 fAti.clear(); 83 94 fDi.clear(); 84 95 fWi.clear(); … … 88 99 { 89 100 fxCst.Clear(); 90 fCovCst.Clear(); 101 fCovCst.Clear(); 102 fCovCstInv.Clear(); 91 103 fXv.Clear(); 92 104 fcovXv.Clear(); … … 95 107 for (Int_t i = 0; i < fNtr; i++) 96 108 { 97 fPar[i]->Clear(); delete fPar[i]; 98 fCov[i]->Clear(); delete fCov[i]; 109 fPar[i]->Clear(); delete fPar[i]; 110 fParNew[i]->Clear(); delete fParNew[i]; 111 fCov[i]->Clear(); delete fCov[i]; 112 fCovNew[i]->Clear(); delete fCovNew[i]; 99 113 } 100 114 fPar.clear(); 101 fCov.clear(); 115 fParNew.clear(); 116 fCov.clear(); 117 fCovNew.clear(); 102 118 // 103 119 ResetWrkArrays(); … … 424 440 // 425 441 // Get track parameters and their covariance 426 TVectorD par = *fPar [i];442 TVectorD par = *fParNew[i]; 427 443 TMatrixDSym Cov = *fCov[i]; 428 444 // … … 433 449 // 434 450 TMatrixD A = Fill_A(par, fs); // A = dx/da = derivatives wrt track parameters 451 TMatrixD At(TMatrixD::kTransposed, A); // A transposed 452 fAti.push_back(new TMatrixD(At)); // Store A' 453 TVectorD di = A * (par - *fPar[i]); // x-shift 454 fdi.push_back(new TVectorD(di)); // Store x-shift 435 455 TMatrixDSym Winv = Cov.Similarity(A); // W^-1 = A*C*A' 436 fWinvi.push_back(new TMatrixDSym(Winv)); 456 fWinvi.push_back(new TMatrixDSym(Winv)); // Store W^-1 matrix 437 457 // 438 458 TMatrixDSym W = RegInv(Winv); // W = (A*C*A')^-1 … … 456 476 { 457 477 //std::cout << "VertexFitter: just in" << std::endl; 458 if (fNtr < 2) 459 { 478 if (fNtr < 2 && !fVtxCst){ 460 479 std::cout << "VertexFit::VertexFitter - Method called with less than 2 tracks - Aborting " << std::endl; 461 480 std::exit(1); … … 469 488 Double_t Chi2 = 0; 470 489 TVectorD x0 = fXv; // If previous fit done 471 if (fRold < 0.0)for (Int_t i = 0; i < 3; i++)x0(i) = 1000.; // Set to arbitrary large value if not 490 if (fRold < 0.0) { 491 // External constraint 492 if (fVtxCst) fRold = TMath::Sqrt(fxCst(0) * fxCst(0) + fxCst(1) * fxCst(1)); 493 // No constraint 494 else for (Int_t i = 0; i < 3; i++)x0(i) = 1000.; // Set to arbitrary large value 495 } 472 496 // 473 497 // Starting vertex radius approximation … … 475 499 Double_t R = fRold; // Use previous fit if available 476 500 if (R < 0.0) R = StartRadius(); // Rough vertex estimate 501 //std::cout << "Start radius: " << R << std::endl; 477 502 // 478 503 // Iteration properties … … 529 554 H += Ds; 530 555 // update constant term 531 TVectorD xs = *fx0i[i]; 556 TVectorD xs = *fx0i[i] - *fdi[i]; 557 TVectorD xx0 = *fx0i[i]; 558 /* 559 std::cout << "Iter. " << Ntry << ", trk " << i << ", x= " 560 << xx0(0) << ", " << xx0(1) << ", " << xx0(2)<< 561 ", ph0= "<<par(1)<< std::endl; 562 */ 532 563 cterm += Ds * xs; 533 564 } // End loop on tracks 565 // Some additions in case of external constraints 566 if (fVtxCst) { 567 H += fCovCstInv; 568 cterm += fCovCstInv * fxCst; 569 DW1D += fCovCstInv; 570 } 534 571 // 535 572 // update vertex position … … 544 581 for (Int_t i = 0; i < fNtr; i++) 545 582 { 546 TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x );583 TVectorD lambda = (*fDi[i]) * (*fx0i[i] - x - *fdi[i]); 547 584 TMatrixDSym Wm1 = *fWinvi[i]; 548 585 fChi2List(i) = Wm1.Similarity(lambda); … … 551 588 TVectorD b = (*fWi[i]) * (x - (*fx0i[i])); 552 589 for (Int_t j = 0; j < 3; j++)ffi[i] += a(j) * b(j) / fa2i[i]; 590 TVectorD newPar = *fPar[i] - ((*fCov[i]) * (*fAti[i])) * lambda; 591 fParNew[i] = new TVectorD(newPar); 553 592 } 593 // Add external constraint to Chi2 594 if (fVtxCst) Chi2 += fCovCstInv.Similarity(x - fxCst); 554 595 // 555 596 TVectorD dx = x - x0; … … 562 603 // Store result 563 604 // 564 fXv = x; 605 fXv = x; // Vertex position 565 606 fcovXv = covX; // Vertex covariance 566 607 fChi2 = Chi2; // Vertex fit Chi2 608 //std::cout << "Found vertex " << fXv(0) << ", " << fXv(1) << ", " << fXv(2) << std::endl; 567 609 } // end of iteration loop 568 610 // … … 603 645 void VertexFit::AddVtxConstraint(TVectorD xv, TMatrixDSym cov) // Add gaussian vertex constraint 604 646 { 605 std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl; 647 //std::cout << "VertexFit::AddVtxConstraint: Not implemented yet" << std::endl; 648 fVtxCst = kTRUE; // Vertex constraint flag 649 fxCst = xv; // Constraint value 650 fCovCst = cov; // Constraint covariance 651 fCovCstInv = cov; 652 fCovCstInv.Invert(); // Its inverse 653 // 654 // Set starting vertex as external constraint 655 fXv = fxCst; 656 fcovXv = fCovCst; 606 657 } 607 658 // … … 613 664 fPar.push_back(par); // add new track 614 665 fCov.push_back(Cov); 666 fParNew.push_back(par); // add new track 667 fCovNew.push_back(Cov); 615 668 // 616 669 // Reset previous vertex temp arrays … … 628 681 fPar.erase(fPar.begin() + iTrk); // Remove track 629 682 fCov.erase(fCov.begin() + iTrk); 683 fParNew.erase(fParNew.begin() + iTrk); // Remove track 684 fCovNew.erase(fCovNew.begin() + iTrk); 630 685 // 631 686 // Reset previous vertex temp arrays -
external/TrackCovariance/VertexFit.h
r952bbbc r53f4746 23 23 Int_t fNtr; // Number of tracks 24 24 std::vector<TVectorD*> fPar; // Input parameter array 25 std::vector<TMatrixDSym*> fCov;// Input parameter covariances 25 std::vector<TVectorD*> fParNew; // Updated parameter array 26 std::vector<TMatrixDSym*> fCov; // Input parameter covariances 27 std::vector<TMatrixDSym*> fCovNew; // Updated parameter covariances 26 28 // Constraints 27 29 Bool_t fVtxCst; // Vertex constraint flag 28 30 TVectorD fxCst; // Constraint value 29 TMatrixDSym fCovCst; // Constraint covariance 31 TMatrixDSym fCovCst; // Constraint 32 TMatrixDSym fCovCstInv; // Inverse of constraint covariance 30 33 // 31 34 // Results … … 39 42 // Work arrays 40 43 std::vector<Double_t> ffi; // Fit phases 41 std::vector<TVectorD*> fx0i; 44 std::vector<TVectorD*> fx0i; // Track expansion points 42 45 std::vector<TVectorD*> fai; // dx/dphi 46 std::vector<TVectorD*> fdi; // x-shift 43 47 std::vector<Double_t> fa2i; // a'Wa 48 std::vector<TMatrixD*> fAti; // A transposed 44 49 std::vector<TMatrixDSym*> fDi; // W-WBW 45 50 std::vector<TMatrixDSym*> fWi; // (ACA')^-1
Note:
See TracChangeset
for help on using the changeset viewer.
