Changeset ebf40fd in git for external/TrackCovariance/ObsTrk.cc
- Timestamp:
- Nov 29, 2021, 3:18:22 PM (3 years ago)
- Branches:
- master
- Children:
- bd376e3
- Parents:
- 9a7ea36
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/ObsTrk.cc
r9a7ea36 rebf40fd 6 6 #include <TRandom.h> 7 7 #include <iostream> 8 #include "SolGeom.h" 8 9 #include "SolGridCov.h" 9 10 #include "ObsTrk.h" … … 11 12 // Constructors 12 13 // 13 // x(3) track origin, p(3) track momentum at origin, Q charge, B ma gnetic field in Tesla14 ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC)14 // x(3) track origin, p(3) track momentum at origin, Q charge, B ma gnetic field in Tesla 15 ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, SolGridCov *GC, SolGeom *G) 15 16 { 16 SetB(B); 17 fB = G->B(); 18 SetB(fB); 19 fG = G; 17 20 fGC = GC; 18 21 fGenX = x; 19 22 fGenP = p; 20 23 fGenQ = Q; 21 fB = B;22 24 fGenPar.ResizeTo(5); 23 25 fGenParMm.ResizeTo(5); … … 36 38 fGenParACTS = ParToACTS(fGenPar); 37 39 fGenParILC = ParToILC(fGenPar); 38 /* 39 cout << "ObsTrk::ObsTrk: fGenPar"; 40 for (Int_t i = 0; i < 5; i++)cout << fGenPar(i) << ", "; 41 cout << endl; 42 */ 43 fObsPar = GenToObsPar(fGenPar, fGC); 40 // 41 fObsPar = GenToObsPar(fGenPar); 44 42 fObsParMm = ParToMm(fObsPar); 45 43 fObsParACTS = ParToACTS(fObsPar); … … 54 52 // 55 53 // x[3] track origin, p[3] track momentum at origin, Q charge, B magnetic field in Tesla 56 ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC)54 ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, SolGridCov* GC, SolGeom *G) 57 55 { 58 SetB(B); 56 fB = G->B(); 57 SetB(fB); 58 fG = G; 59 59 fGC = GC; 60 60 fGenX.SetXYZ(x[0],x[1],x[2]); 61 61 fGenP.SetXYZ(p[0],p[1],p[2]); 62 62 fGenQ = Q; 63 fB = B;64 63 fGenPar.ResizeTo(5); 65 64 fGenParMm.ResizeTo(5); … … 77 76 fGenParACTS = ParToACTS(fGenPar); 78 77 fGenParILC = ParToILC(fGenPar); 79 /* 80 cout << "ObsTrk::ObsTrk: fGenPar"; 81 for (Int_t i = 0; i < 5; i++)cout << fGenPar(i) << ", "; 82 cout << endl; 83 */ 84 fObsPar = GenToObsPar(fGenPar, fGC); 78 // 79 fObsPar = GenToObsPar(fGenPar); 85 80 fObsParACTS = ParToACTS(fObsPar); 86 81 fObsParILC = ParToILC(fObsPar); … … 113 108 } 114 109 // 115 TVectorD ObsTrk::GenToObsPar(TVectorD gPar , SolGridCov *GC)110 TVectorD ObsTrk::GenToObsPar(TVectorD gPar) 116 111 { 117 TVector3 p = ParToP(gPar);118 Double_t pt = p.Pt();119 Double_t tanTh = 1.0 / TMath::Abs(gPar(4));120 Double_t angd = TMath::ATan(tanTh)*180. / TMath::Pi();121 112 // 122 113 // Check ranges 123 Double_t minPt = GC->GetMinPt();124 //if (pt < minPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt <<endl;125 Double_t maxPt = GC->GetMaxPt();126 //if (pt > maxPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt <<endl;127 Double_t minAn = GC->GetMinAng();128 //if (angd < minAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd129 // << " is below grid range of " << minAn << endl;130 Double_t maxAn = GC->GetMaxAng();131 //if (angd > maxAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd132 // << " is above grid range of " << maxAn << endl;114 Double_t minPt = fGC->GetMinPt(); 115 //if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl; 116 Double_t maxPt = fGC->GetMaxPt(); 117 //if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl; 118 Double_t minAn = fGC->GetMinAng(); 119 //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 120 // << " is below grid range of " << minAn << std::endl; 121 Double_t maxAn = fGC->GetMaxAng(); 122 //if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 123 // << " is above grid range of " << maxAn << std::endl; 133 124 // 134 TMatrixDSym Cov = GC->GetCov(pt, angd); 125 TMatrixDSym Cov(5); 126 // 127 // Check if track origin is inside beampipe and betwen the first disks 128 // 129 Double_t Rin = fG->GetRmin(); 130 Double_t ZinPos = fG->GetZminPos(); 131 Double_t ZinNeg = fG->GetZminNeg(); 132 Bool_t inside = TrkUtil::IsInside(fGenX, Rin, ZinNeg, ZinPos); // Check if in inner box 133 if (inside) 134 { 135 //std::cout<<"ObsTrk:: inside: x= "<<fGenX(0)<<", y= "<<fGenX(1) 136 // <<", z= "<<fGenX(2)<<std::endl; 137 // Observed track parameters 138 Double_t pt = fGenP.Pt(); 139 Double_t angd = fGenP.Theta() * 180. / TMath::Pi(); 140 Cov = fGC->GetCov(pt, angd); // Track covariance 141 } 142 else 143 { 144 //std::cout<<"ObsTrk:: outside: x= "<<fGenX(0)<<", y= "<<fGenX(1) 145 // <<", z= "<<fGenX(2)<<std::endl; 146 SolTrack* trk = new SolTrack(fGenX, fGenP, fG); 147 Bool_t Res = kTRUE; Bool_t MS = kTRUE; 148 trk->CovCalc(Res, MS); // Calculate covariance matrix 149 Cov = trk->Cov(); // Track covariance 150 delete trk; 151 } 152 // 135 153 fCov = Cov; 136 154 // 137 155 // Now do Choleski decomposition and random number extraction, with appropriate stabilization 138 156 // 139 TMatrixDSym CvN = Cov; 140 TMatrixDSym DCv(5); DCv.Zero(); 141 TMatrixDSym DCvInv(5); DCvInv.Zero(); 142 for (Int_t id = 0; id < 5; id++) 143 { 144 Double_t dVal = TMath::Sqrt(Cov(id, id)); 145 DCv (id, id) = dVal; 146 DCvInv(id, id) = 1.0 / dVal; 147 } 148 CvN.Similarity(DCvInv); // Normalize diagonal to 1 149 TDecompChol Chl(CvN); 150 Bool_t OK = Chl.Decompose(); // Choleski decomposition of normalized matrix 151 TMatrixD U = Chl.GetU(); // Get Upper triangular matrix 152 TMatrixD Ut(TMatrixD::kTransposed, U); // Transposed of U (lower triangular) 153 TVectorD r(5); 154 for (Int_t i = 0; i < 5; i++)r(i) = gRandom->Gaus(0.0, 1.0); // Array of normal random numbers 155 TVectorD oPar = gPar + DCv*(Ut*r); // Observed parameter vector 157 TVectorD oPar = TrkUtil::CovSmear(gPar, Cov); 156 158 // 157 159 return oPar;
Note:
See TracChangeset
for help on using the changeset viewer.