Changeset a617744 in git
- Timestamp:
- Mar 1, 2021, 4:01:37 PM (4 years ago)
- Branches:
- master
- Children:
- 40e890c
- Parents:
- 3f8466a
- Location:
- external/TrackCovariance
- Files:
-
- 2 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
external/TrackCovariance/ObsTrk.cc
r3f8466a 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; 36 42 */ 37 43 fObsPar = GenToObsPar(fGenPar, fGC); 44 fObsParMm = ParToMm(fObsPar); 38 45 fObsParACTS = ParToACTS(fObsPar); 39 46 fObsParILC = ParToILC(fObsPar); … … 41 48 fObsP = ParToP(fObsPar); 42 49 fObsQ = ParToQ(fObsPar); 43 fCovACTS = CovToACTS(fCov); 50 fCovMm = CovToMm(fCov); 51 fCovACTS = CovToACTS(fObsPar, fCov); 44 52 fCovILC = CovToILC(fCov); 45 53 } 46 54 // 47 55 // x[3] track origin, p[3] track momentum at origin, Q charge, B magnetic field in Tesla 48 56 ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC) 49 57 { 58 SetBfield(B); 50 59 fGC = GC; 51 60 fGenX.SetXYZ(x[0],x[1],x[2]); … … 54 63 fB = B; 55 64 fGenPar.ResizeTo(5); 65 fGenParMm.ResizeTo(5); 56 66 fGenParACTS.ResizeTo(6); 57 67 fGenParILC.ResizeTo(5); 58 68 fObsPar.ResizeTo(5); 69 fObsParMm.ResizeTo(5); 59 70 fObsParACTS.ResizeTo(6); 60 71 fObsParILC.ResizeTo(5); 61 72 fCov.ResizeTo(5, 5); 73 fCovMm.ResizeTo(5, 5); 62 74 fCovACTS.ResizeTo(6, 6); 63 75 fCovILC.ResizeTo(5, 5); … … 76 88 fObsP = ParToP(fObsPar); 77 89 fObsQ = ParToQ(fObsPar); 78 fCovACTS = CovToACTS(f Cov);90 fCovACTS = CovToACTS(fObsPar, fCov); 79 91 fCovILC = CovToILC(fCov); 80 92 } 81 82 83 93 // 84 94 // Destructor … … 88 98 fGenP.Clear(); 89 99 fGenPar.Clear(); 100 fGenParMm.Clear(); 90 101 fGenParACTS.Clear(); 102 fGenParILC.Clear(); 91 103 fObsX.Clear(); 92 104 fObsP.Clear(); 93 105 fObsPar.Clear(); 106 fObsParMm.Clear(); 94 107 fObsParACTS.Clear(); 108 fObsParILC.Clear(); 95 109 fCov.Clear(); 110 fCovMm.Clear(); 96 111 fCovACTS.Clear(); 97 } 98 TVectorD ObsTrk::XPtoPar(TVector3 x, TVector3 p, Double_t Q) 99 { 100 // 101 TVectorD Par(5); 102 // Transverse parameters 103 Double_t a = -Q*fB*0.2998; // Units are Tesla, GeV and meters 104 Double_t pt = p.Pt(); 105 Double_t C = a / (2 * pt); // Half curvature 106 //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl; 107 Double_t r2 = x.Perp2(); 108 Double_t cross = x(0)*p(1) - x(1)*p(0); 109 Double_t T = TMath::Sqrt(pt*pt - 2 * a*cross + a*a*r2); 110 Double_t phi0 = TMath::ATan2((p(1) - a*x(0)) / T, (p(0) + a*x(1)) / T); // Phi0 111 Double_t D; // Impact parameter D 112 if (pt < 10.0) D = (T - pt) / a; 113 else D = (-2 * cross + a*r2) / (T + pt); 114 // 115 Par(0) = D; // Store D 116 Par(1) = phi0; // Store phi0 117 Par(2) = C; // Store C 118 //Longitudinal parameters 119 Double_t B = C*TMath::Sqrt(TMath::Max(r2 - D*D,0.0) / (1 + 2 * C*D)); 120 Double_t st = TMath::ASin(B) / C; 121 Double_t ct = p(2) / pt; 122 Double_t z0 = x(2) - ct*st; 123 // 124 Par(3) = z0; // Store z0 125 Par(4) = ct; // Store cot(theta) 126 // 127 return Par; 128 } 129 // 130 TVector3 ObsTrk::ParToX(TVectorD Par) 131 { 132 Double_t D = Par(0); 133 Double_t phi0 = Par(1); 134 Double_t z0 = Par(3); 135 // 136 TVector3 Xval; 137 Xval(0) = -D*TMath::Sin(phi0); 138 Xval(1) = D*TMath::Cos(phi0); 139 Xval(2) = z0; 140 // 141 return Xval; 142 } 143 // 144 TVector3 ObsTrk::ParToP(TVectorD Par) 145 { 146 Double_t C = Par(2); 147 Double_t phi0 = Par(1); 148 Double_t ct = Par(4); 149 // 150 TVector3 Pval; 151 Double_t pt = fB*0.2998 / TMath::Abs(2 * C); 152 Pval(0) = pt*TMath::Cos(phi0); 153 Pval(1) = pt*TMath::Sin(phi0); 154 Pval(2) = pt*ct; 155 // 156 return Pval; 157 } 158 // 159 160 Double_t ObsTrk::ParToQ(TVectorD Par) 161 { 162 return TMath::Sign(1.0, -Par(2)); 112 fCovILC.Clear(); 163 113 } 164 114 // … … 172 122 // Check ranges 173 123 Double_t minPt = GC->GetMinPt (); 174 //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; 175 125 Double_t maxPt = GC->GetMaxPt(); 176 //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; 177 127 Double_t minAn = GC->GetMinAng(); 178 //if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 179 // << " 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; 180 130 Double_t maxAn = GC->GetMaxAng(); 181 //if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd182 // << " 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; 183 133 // 184 134 TMatrixDSym Cov = GC->GetCov(pt, angd); … … 207 157 return oPar; 208 158 } 209 // Parameter conversion to ACTS format210 TVectorD ObsTrk::ParToACTS(TVectorD Par)211 {212 TVectorD pACTS(6); // Return vector213 //214 Double_t b = -0.29988*fB / 2.;215 pACTS(0) = 1000*Par(0); // D from m to mm216 pACTS(1) = 1000 * Par(3); // z0 from m to mm217 pACTS(2) = Par(1); // Phi0 is unchanged218 pACTS(3) = TMath::ATan2(1.0,Par(4)); // Theta in [0, pi] range219 pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4))); // q/p in GeV220 pACTS(5) = 0.0; // Time: currently undefined221 //222 return pACTS;223 }224 // Covariance conversion to ACTS format225 TMatrixDSym ObsTrk::CovToACTS(TMatrixDSym Cov)226 {227 TMatrixDSym cACTS(6); cACTS.Zero();228 Double_t b = -0.29988*fB / 2.;229 //230 // Fill derivative matrix231 TMatrixD A(5, 5); A.Zero();232 Double_t ct = fGenPar(4); // cot(theta)233 Double_t C = fGenPar(2); // half curvature234 A(0, 0) = 1000.; // D-D conversion to mm235 A(1, 2) = 1.0; // phi0-phi0236 A(2, 4) = 1.0/(TMath::Sqrt(1.0 + ct*ct) * b); // q/p-C237 A(3, 1) = 1000.; // z0-z0 conversion to mm238 A(4, 3) = -1.0 / (1.0 + ct*ct); // theta - cot(theta)239 A(4, 4) = -C*ct / (b*pow(1.0 + ct*ct,3.0/2.0)); // q/p-cot(theta)240 //241 TMatrixDSym Cv = Cov;242 TMatrixD At(5, 5);243 At.Transpose(A);244 Cv.Similarity(At);245 TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv;246 cACTS(5, 5) = 0.1; // Currently undefined: set to arbitrary value to avoid crashes247 //248 return cACTS;249 }250 251 // Parameter conversion to ILC format252 TVectorD ObsTrk::ParToILC(TVectorD Par)253 {254 TVectorD pILC(5); // Return vector255 //256 pILC(0) = Par(0)*1.0e3; // d0 in mm257 pILC(1) = Par(1); // phi0 is unchanged258 pILC(2) = -2 * Par(2)*1.0e-3; // w in mm^-1259 pILC(3) = Par(3)*1.0e3; // z0 in mm260 pILC(4) = Par(4); // tan(lambda) = cot(theta)261 //262 return pILC;263 }264 // Covariance conversion to ILC format265 TMatrixDSym ObsTrk::CovToILC(TMatrixDSym Cov)266 {267 TMatrixDSym cILC(5); cILC.Zero();268 //269 // Fill derivative matrix270 TMatrixD A(5, 5); A.Zero();271 //272 A(0, 0) = 1.0e3; // D-d0 in mm273 A(1, 1) = 1.0; // phi0-phi0274 A(2, 2) = -2.0e-3; // w-C275 A(3, 3) = 1.0e3; // z0-z0 conversion to mm276 A(4, 4) = 1.0; // tan(lambda) - cot(theta)277 //278 TMatrixDSym Cv = Cov;279 TMatrixD At(5, 5);280 At.Transpose(A);281 Cv.Similarity(At);282 cILC = Cv;283 //284 return cILC;285 } -
external/TrackCovariance/ObsTrk.h
r3f8466a ra617744 6 6 #include <TMatrixDSym.h> 7 7 #include <TDecompChol.h> 8 #include "TrkUtil.h" 8 9 #include "SolGridCov.h" 9 10 // … … 13 14 // INFN - Sezione di Pisa, Italy 14 15 // 15 class ObsTrk{ 16 class ObsTrk: public TrkUtil 17 { 16 18 // 17 19 // Class to handle simulation of tracking resolution … … 19 21 // Prefix Gen marks variables before resolution smearing 20 22 // 21 private: 22 Double_t fB; 23 SolGridCov *fGC; 23 private: 24 Double_t fB; // Solenoid magnetic field 25 SolGridCov *fGC; // Covariance matrix grid 24 26 Double_t fGenQ; // Generated track charge 25 27 Double_t fObsQ; // Observed track charge 26 28 TVector3 fGenX; // Generated track origin (x,y,z) 27 TVector3 fObsX; // Observed track origin (x,y,z) @ track min. approach 29 TVector3 fObsX; // Observed track origin (x,y,z) @ track min. approach 28 30 TVector3 fGenP; // Generated track momentum at track origin 29 31 TVector3 fObsP; // Observed track momentum @ track minimum approach 30 TVectorD fGenPar; // Generated helix track parameters (D, phi0, C, z0, cot(th)) 32 TVectorD fGenPar; // Generated helix track parameters (D, phi0, C, z0, cot(th)) in meters 33 TVectorD fGenParMm; // Generated helix track parameters (D, phi0, C, z0, cot(th)) in mm 31 34 TVectorD fGenParACTS; // Generated helix track parameters (D, z0, phi0, th, q/p, time 32 TVectorD fGenParILC; // Generated helix track parameters (w, phi0, d0, z0, tan(lambda)) 33 TVectorD fObsPar; // Observed helix track parameters (D, phi0, C, z0, cot(th)) 35 TVectorD fGenParILC; // Generated helix track parameters (w, phi0, d0, z0, tan(lambda)) 36 TVectorD fObsPar; // Observed helix track parameters (D, phi0, C, z0, cot(th)) in meters 37 TVectorD fObsParMm; // Observed helix track parameters (D, phi0, C, z0, cot(th)) in mm 34 38 TVectorD fObsParACTS; // Observed helix track parameters (D, z0, phi0, th, q/p, time 35 TVectorD fObsParILC; // Observed helix track parameters (d0, phi0, w, z0, tan(lambda)) 36 TMatrixDSym fCov; // INterpolated covariance of track parameters 39 TVectorD fObsParILC; // Observed helix track parameters (d0, phi0, w, z0, tan(lambda)) 40 TMatrixDSym fCov; // Interpolated covariance of track in meters 41 TMatrixDSym fCovMm; // Interpolated covariance of track parameters in mm 37 42 TMatrixDSym fCovACTS; // Covariance of track parameters in ACTS format 38 43 // (D, z0, phi0, theta, q/p, time) 39 TMatrixDSym fCovILC; 44 TMatrixDSym fCovILC; // Covariance of track parameters in ILC format 40 45 // (d0, phi0, w, z0, tan(lambda)) 41 46 // 42 // Conversion to ACTS parametrization47 // Service routines 43 48 // 44 TVectorD ParToACTS(TVectorD Par); // Parameter conversion 45 TMatrixDSym CovToACTS(TMatrixDSym Cov); // Covariance 46 // 47 // Conversion to ILC parametrization 48 // 49 TVectorD ParToILC(TVectorD Par); // Parameter conversion 50 TMatrixDSym CovToILC(TMatrixDSym Cov); // Covariance conversion 49 TVectorD GenToObsPar(TVectorD gPar, SolGridCov* GC); 51 50 // 52 51 public: … … 54 53 // Constructors 55 54 // x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla 56 ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC); // Initialize and generate smeared track55 ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC); // Initialize and generate smeared 57 56 ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC); // Initialize and generate smeared track 58 57 // Destructor 59 58 ~ObsTrk(); 60 //61 // Service routines62 //63 TVectorD XPtoPar(TVector3 x, TVector3 p, Double_t Q);64 TVectorD GenToObsPar(TVectorD gPar, SolGridCov *GC);65 TVector3 ParToX(TVectorD Par);66 TVector3 ParToP(TVectorD Par);67 Double_t ParToQ(TVectorD Par);68 59 // 69 60 // Accessors … … 75 66 TVector3 GetGenP() { return fGenP; } 76 67 // D, phi0, C, z0, cot(th) 77 TVectorD GetGenPar() { return fGenPar; } 68 TVectorD GetGenPar() { return fGenPar; } // in meters 69 TVectorD GetGenParMm() { return fGenParMm; } // in mm 78 70 // D, z0, phi0, theta, q/p, time 79 71 TVectorD GetGenParACTS() { return fGenParACTS; } … … 85 77 TVector3 GetObsP() { return fObsP; } 86 78 // D, phi0, C, z0, cot(th) 87 TVectorD GetObsPar() { return fObsPar; } 79 TVectorD GetObsPar() { return fObsPar; } // in meters 80 TVectorD GetObsParMm() { return fObsParMm; } // In mm 88 81 // D, z0, phi0, theta, q/p, time 89 82 TVectorD GetObsParACTS() { return fObsParACTS; } … … 91 84 TVectorD GetObsParILC() { return fObsParILC; } 92 85 // Covariances 93 TMatrixDSym GetCov(){ return fCov; } 86 TMatrixDSym GetCov() { return fCov; } // in meters 87 TMatrixDSym GetCovMm() { return fCov; } // in mm 94 88 TMatrixDSym GetCovACTS(){ return fCovACTS; } 95 89 TMatrixDSym GetCovILC(){ return fCovILC; } 90 // 96 91 }; 97 92
Note:
See TracChangeset
for help on using the changeset viewer.