Changes in external/TrackCovariance/ObsTrk.cc [a0f5d71:a617744] in git
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/ObsTrk.cc
ra0f5d71 ra617744 10 10 // 11 11 // Constructors 12 // 12 13 // x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla 13 14 ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC) 14 15 { 16 SetBfield(B); 15 17 fGC = GC; 16 18 fGenX = x; … … 19 21 fB = B; 20 22 fGenPar.ResizeTo(5); 23 fGenParMm.ResizeTo(5); 21 24 fGenParACTS.ResizeTo(6); 22 25 fGenParILC.ResizeTo(5); 23 26 fObsPar.ResizeTo(5); 27 fObsParMm.ResizeTo(5); 24 28 fObsParACTS.ResizeTo(6); 25 29 fObsParILC.ResizeTo(5); 26 30 fCov.ResizeTo(5, 5); 31 fCovMm.ResizeTo(5, 5); 27 32 fCovACTS.ResizeTo(6, 6); 28 33 fCovILC.ResizeTo(5, 5); 29 34 fGenPar = XPtoPar(x,p,Q); 35 fGenParMm = ParToMm(fGenPar); 30 36 fGenParACTS = ParToACTS(fGenPar); 31 37 fGenParILC = ParToILC(fGenPar); 32 38 /* 33 std::cout << "ObsTrk::ObsTrk: fGenPar"; 34 for (Int_t i = 0; i < 5; i++)std::cout << fGenPar(i) << ", "; 35 std::cout << std::endl; 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); 44 fObsParMm = ParToMm(fObsPar); 45 fObsParACTS = ParToACTS(fObsPar); 46 fObsParILC = ParToILC(fObsPar); 47 fObsX = ParToX(fObsPar); 48 fObsP = ParToP(fObsPar); 49 fObsQ = ParToQ(fObsPar); 50 fCovMm = CovToMm(fCov); 51 fCovACTS = CovToACTS(fObsPar, fCov); 52 fCovILC = CovToILC(fCov); 53 } 54 // 55 // 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) 57 { 58 SetBfield(B); 59 fGC = GC; 60 fGenX.SetXYZ(x[0],x[1],x[2]); 61 fGenP.SetXYZ(p[0],p[1],p[2]); 62 fGenQ = Q; 63 fB = B; 64 fGenPar.ResizeTo(5); 65 fGenParMm.ResizeTo(5); 66 fGenParACTS.ResizeTo(6); 67 fGenParILC.ResizeTo(5); 68 fObsPar.ResizeTo(5); 69 fObsParMm.ResizeTo(5); 70 fObsParACTS.ResizeTo(6); 71 fObsParILC.ResizeTo(5); 72 fCov.ResizeTo(5, 5); 73 fCovMm.ResizeTo(5, 5); 74 fCovACTS.ResizeTo(6, 6); 75 fCovILC.ResizeTo(5, 5); 76 fGenPar = XPtoPar(fGenX, fGenP, Q); 77 fGenParACTS = ParToACTS(fGenPar); 78 fGenParILC = ParToILC(fGenPar); 79 /* 80 cout << "ObsTrk::ObsTrk: fGenPar"; 81 for (Int_t i = 0; i < 5; i++)cout << fGenPar(i) << ", "; 82 cout << endl; 36 83 */ 37 84 fObsPar = GenToObsPar(fGenPar, fGC); … … 41 88 fObsP = ParToP(fObsPar); 42 89 fObsQ = ParToQ(fObsPar); 43 fCovACTS = CovToACTS(f Cov);90 fCovACTS = CovToACTS(fObsPar, fCov); 44 91 fCovILC = CovToILC(fCov); 45 92 } … … 51 98 fGenP.Clear(); 52 99 fGenPar.Clear(); 100 fGenParMm.Clear(); 53 101 fGenParACTS.Clear(); 102 fGenParILC.Clear(); 54 103 fObsX.Clear(); 55 104 fObsP.Clear(); 56 105 fObsPar.Clear(); 106 fObsParMm.Clear(); 57 107 fObsParACTS.Clear(); 108 fObsParILC.Clear(); 58 109 fCov.Clear(); 110 fCovMm.Clear(); 59 111 fCovACTS.Clear(); 60 } 61 TVectorD ObsTrk::XPtoPar(TVector3 x, TVector3 p, Double_t Q) 62 { 63 // 64 TVectorD Par(5); 65 // Transverse parameters 66 Double_t a = -Q*fB*0.2998; // Units are Tesla, GeV and meters 67 Double_t pt = p.Pt(); 68 Double_t C = a / (2 * pt); // Half curvature 69 //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl; 70 Double_t r2 = x.Perp2(); 71 Double_t cross = x(0)*p(1) - x(1)*p(0); 72 Double_t T = TMath::Sqrt(pt*pt - 2 * a*cross + a*a*r2); 73 Double_t phi0 = TMath::ATan2((p(1) - a*x(0)) / T, (p(0) + a*x(1)) / T); // Phi0 74 Double_t D; // Impact parameter D 75 if (pt < 10.0) D = (T - pt) / a; 76 else D = (-2 * cross + a*r2) / (T + pt); 77 // 78 Par(0) = D; // Store D 79 Par(1) = phi0; // Store phi0 80 Par(2) = C; // Store C 81 //Longitudinal parameters 82 Double_t B = C*TMath::Sqrt(TMath::Max(r2 - D*D,0.0) / (1 + 2 * C*D)); 83 Double_t st = TMath::ASin(B) / C; 84 Double_t ct = p(2) / pt; 85 Double_t z0 = x(2) - ct*st; 86 // 87 Par(3) = z0; // Store z0 88 Par(4) = ct; // Store cot(theta) 89 // 90 return Par; 91 } 92 // 93 TVector3 ObsTrk::ParToX(TVectorD Par) 94 { 95 Double_t D = Par(0); 96 Double_t phi0 = Par(1); 97 Double_t z0 = Par(3); 98 // 99 TVector3 Xval; 100 Xval(0) = -D*TMath::Sin(phi0); 101 Xval(1) = D*TMath::Cos(phi0); 102 Xval(2) = z0; 103 // 104 return Xval; 105 } 106 // 107 TVector3 ObsTrk::ParToP(TVectorD Par) 108 { 109 Double_t C = Par(2); 110 Double_t phi0 = Par(1); 111 Double_t ct = Par(4); 112 // 113 TVector3 Pval; 114 Double_t pt = fB*0.2998 / TMath::Abs(2 * C); 115 Pval(0) = pt*TMath::Cos(phi0); 116 Pval(1) = pt*TMath::Sin(phi0); 117 Pval(2) = pt*ct; 118 // 119 return Pval; 120 } 121 // 122 123 Double_t ObsTrk::ParToQ(TVectorD Par) 124 { 125 return TMath::Sign(1.0, -Par(2)); 112 fCovILC.Clear(); 126 113 } 127 114 // … … 135 122 // Check ranges 136 123 Double_t minPt = GC->GetMinPt (); 137 if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;124 //if (pt < minPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << endl; 138 125 Double_t maxPt = GC->GetMaxPt(); 139 if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;126 //if (pt > maxPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << endl; 140 127 Double_t minAn = GC->GetMinAng(); 141 if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd142 << " is below grid range of " << minAn << std::endl;128 //if (angd < minAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd 129 // << " is below grid range of " << minAn << endl; 143 130 Double_t maxAn = GC->GetMaxAng(); 144 if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd145 << " is above grid range of " << maxAn << std::endl;131 //if (angd > maxAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd 132 // << " is above grid range of " << maxAn << endl; 146 133 // 147 134 TMatrixDSym Cov = GC->GetCov(pt, angd); … … 170 157 return oPar; 171 158 } 172 // Parameter conversion to ACTS format173 TVectorD ObsTrk::ParToACTS(TVectorD Par)174 {175 TVectorD pACTS(6); // Return vector176 //177 Double_t b = -0.29988*fB / 2.;178 pACTS(0) = 1000*Par(0); // D from m to mm179 pACTS(1) = 1000 * Par(3); // z0 from m to mm180 pACTS(2) = Par(1); // Phi0 is unchanged181 pACTS(3) = TMath::ATan(1.0 / Par(4)) + TMath::PiOver2(); // Theta in [0, pi] range182 pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4))); // q/p in GeV183 pACTS(5) = 0.0; // Time: currently undefined184 //185 return pACTS;186 }187 // Covariance conversion to ACTS format188 TMatrixDSym ObsTrk::CovToACTS(TMatrixDSym Cov)189 {190 TMatrixDSym cACTS(6); cACTS.Zero();191 Double_t b = -0.29988*fB / 2.;192 //193 // Fill derivative matrix194 TMatrixD A(5, 5); A.Zero();195 Double_t ct = fGenPar(4); // cot(theta)196 Double_t C = fGenPar(2); // half curvature197 A(0, 0) = 1000.; // D-D conversion to mm198 A(1, 2) = 1.0; // phi0-phi0199 A(2, 4) = 1.0/(TMath::Sqrt(1.0 + ct*ct) * b); // q/p-C200 A(3, 1) = 1000.; // z0-z0 conversion to mm201 A(4, 3) = -1.0 / (1.0 + ct*ct); // theta - cot(theta)202 A(4, 4) = -C*ct / (b*pow(1.0 + ct*ct,3.0/2.0)); // q/p-cot(theta)203 //204 TMatrixDSym Cv = Cov;205 TMatrixD At(5, 5);206 At.Transpose(A);207 Cv.Similarity(At);208 TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv;209 cACTS(5, 5) = 0.1; // Currently undefined: set to arbitrary value to avoid crashes210 //211 return cACTS;212 }213 214 // Parameter conversion to ILC format215 TVectorD ObsTrk::ParToILC(TVectorD Par)216 {217 TVectorD pILC(5); // Return vector218 //219 pILC(0) = Par(0)*1.0e3; // d0 in mm220 pILC(1) = Par(1); // phi0 is unchanged221 pILC(2) = -2 * Par(2)*1.0e-3; // w in mm^-1222 pILC(3) = Par(3)*1.0e3; // z0 in mm223 pILC(4) = Par(4); // tan(lambda) = cot(theta)224 //225 return pILC;226 }227 // Covariance conversion to ILC format228 TMatrixDSym ObsTrk::CovToILC(TMatrixDSym Cov)229 {230 TMatrixDSym cILC(5); cILC.Zero();231 //232 // Fill derivative matrix233 TMatrixD A(5, 5); A.Zero();234 //235 A(0, 0) = 1.0e3; // D-d0 in mm236 A(1, 1) = 1.0; // phi0-phi0237 A(2, 2) = -2.0e-3; // w-C238 A(3, 3) = 1.0e3; // z0-z0 conversion to mm239 A(4, 4) = 1.0; // tan(lambda) - cot(theta)240 //241 TMatrixDSym Cv = Cov;242 TMatrixD At(5, 5);243 At.Transpose(A);244 Cv.Similarity(At);245 cILC = Cv;246 //247 return cILC;248 }249 250 251 252
Note:
See TracChangeset
for help on using the changeset viewer.