Changes in external/TrackCovariance/ObsTrk.cc [a617744:a0f5d71] in git
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/ObsTrk.cc
ra617744 ra0f5d71 10 10 // 11 11 // Constructors 12 //13 12 // x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla 14 13 ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC) 15 14 { 16 SetBfield(B);17 15 fGC = GC; 18 16 fGenX = x; … … 21 19 fB = B; 22 20 fGenPar.ResizeTo(5); 23 fGenParMm.ResizeTo(5);24 21 fGenParACTS.ResizeTo(6); 25 22 fGenParILC.ResizeTo(5); 26 23 fObsPar.ResizeTo(5); 27 fObsParMm.ResizeTo(5);28 24 fObsParACTS.ResizeTo(6); 29 25 fObsParILC.ResizeTo(5); 30 26 fCov.ResizeTo(5, 5); 31 fCovMm.ResizeTo(5, 5);32 27 fCovACTS.ResizeTo(6, 6); 33 28 fCovILC.ResizeTo(5, 5); 34 29 fGenPar = XPtoPar(x,p,Q); 35 fGenParMm = ParToMm(fGenPar);36 30 fGenParACTS = ParToACTS(fGenPar); 37 31 fGenParILC = ParToILC(fGenPar); 38 32 /* 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; 33 std::cout << "ObsTrk::ObsTrk: fGenPar"; 34 for (Int_t i = 0; i < 5; i++)std::cout << fGenPar(i) << ", "; 35 std::cout << std::endl; 83 36 */ 84 37 fObsPar = GenToObsPar(fGenPar, fGC); … … 88 41 fObsP = ParToP(fObsPar); 89 42 fObsQ = ParToQ(fObsPar); 90 fCovACTS = CovToACTS(f ObsPar, fCov);43 fCovACTS = CovToACTS(fCov); 91 44 fCovILC = CovToILC(fCov); 92 45 } … … 98 51 fGenP.Clear(); 99 52 fGenPar.Clear(); 100 fGenParMm.Clear();101 53 fGenParACTS.Clear(); 102 fGenParILC.Clear();103 54 fObsX.Clear(); 104 55 fObsP.Clear(); 105 56 fObsPar.Clear(); 106 fObsParMm.Clear();107 57 fObsParACTS.Clear(); 108 fObsParILC.Clear();109 58 fCov.Clear(); 110 fCovMm.Clear();111 59 fCovACTS.Clear(); 112 fCovILC.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)); 113 126 } 114 127 // … … 122 135 // Check ranges 123 136 Double_t minPt = GC->GetMinPt (); 124 //if (pt < minPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt <<endl;137 if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl; 125 138 Double_t maxPt = GC->GetMaxPt(); 126 //if (pt > maxPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt <<endl;139 if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl; 127 140 Double_t minAn = GC->GetMinAng(); 128 //if (angd < minAn)cout << "Warning ObsTrk::GenToObsPar: angle " << angd129 // << " is below grid range of " << minAn <<endl;141 if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 142 << " is below grid range of " << minAn << std::endl; 130 143 Double_t maxAn = GC->GetMaxAng(); 131 //if (angd > maxAn)cout << "Warning ObsTrk::GenToObsPar: angle " << angd132 // << " is above grid range of " << maxAn <<endl;144 if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 145 << " is above grid range of " << maxAn << std::endl; 133 146 // 134 147 TMatrixDSym Cov = GC->GetCov(pt, angd); … … 157 170 return oPar; 158 171 } 172 // Parameter conversion to ACTS format 173 TVectorD ObsTrk::ParToACTS(TVectorD Par) 174 { 175 TVectorD pACTS(6); // Return vector 176 // 177 Double_t b = -0.29988*fB / 2.; 178 pACTS(0) = 1000*Par(0); // D from m to mm 179 pACTS(1) = 1000 * Par(3); // z0 from m to mm 180 pACTS(2) = Par(1); // Phi0 is unchanged 181 pACTS(3) = TMath::ATan(1.0 / Par(4)) + TMath::PiOver2(); // Theta in [0, pi] range 182 pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4))); // q/p in GeV 183 pACTS(5) = 0.0; // Time: currently undefined 184 // 185 return pACTS; 186 } 187 // Covariance conversion to ACTS format 188 TMatrixDSym ObsTrk::CovToACTS(TMatrixDSym Cov) 189 { 190 TMatrixDSym cACTS(6); cACTS.Zero(); 191 Double_t b = -0.29988*fB / 2.; 192 // 193 // Fill derivative matrix 194 TMatrixD A(5, 5); A.Zero(); 195 Double_t ct = fGenPar(4); // cot(theta) 196 Double_t C = fGenPar(2); // half curvature 197 A(0, 0) = 1000.; // D-D conversion to mm 198 A(1, 2) = 1.0; // phi0-phi0 199 A(2, 4) = 1.0/(TMath::Sqrt(1.0 + ct*ct) * b); // q/p-C 200 A(3, 1) = 1000.; // z0-z0 conversion to mm 201 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 crashes 210 // 211 return cACTS; 212 } 213 214 // Parameter conversion to ILC format 215 TVectorD ObsTrk::ParToILC(TVectorD Par) 216 { 217 TVectorD pILC(5); // Return vector 218 // 219 pILC(0) = Par(0)*1.0e3; // d0 in mm 220 pILC(1) = Par(1); // phi0 is unchanged 221 pILC(2) = -2 * Par(2)*1.0e-3; // w in mm^-1 222 pILC(3) = Par(3)*1.0e3; // z0 in mm 223 pILC(4) = Par(4); // tan(lambda) = cot(theta) 224 // 225 return pILC; 226 } 227 // Covariance conversion to ILC format 228 TMatrixDSym ObsTrk::CovToILC(TMatrixDSym Cov) 229 { 230 TMatrixDSym cILC(5); cILC.Zero(); 231 // 232 // Fill derivative matrix 233 TMatrixD A(5, 5); A.Zero(); 234 // 235 A(0, 0) = 1.0e3; // D-d0 in mm 236 A(1, 1) = 1.0; // phi0-phi0 237 A(2, 2) = -2.0e-3; // w-C 238 A(3, 3) = 1.0e3; // z0-z0 conversion to mm 239 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.