Changes in external/TrackCovariance/ObsTrk.cc [a617744:a0f5d71] in git

File:

• external/TrackCovariance/ObsTrk.cc (modified) (6 diffs)

external/TrackCovariance/ObsTrk.cc

@@ -10 +10 @@
 //
 // Constructors
-//
 // x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla
 ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC)
 {
-        SetBfield(B);
         fGC = GC;
         fGenX = x;
@@ -21 +19 @@
         fB = B;
         fGenPar.ResizeTo(5);
-        fGenParMm.ResizeTo(5);
         fGenParACTS.ResizeTo(6);
         fGenParILC.ResizeTo(5);
         fObsPar.ResizeTo(5);
-        fObsParMm.ResizeTo(5);
         fObsParACTS.ResizeTo(6);
         fObsParILC.ResizeTo(5);
         fCov.ResizeTo(5, 5);
-        fCovMm.ResizeTo(5, 5);
         fCovACTS.ResizeTo(6, 6);
         fCovILC.ResizeTo(5, 5);
         fGenPar = XPtoPar(x,p,Q);
-        fGenParMm = ParToMm(fGenPar);
         fGenParACTS = ParToACTS(fGenPar);
         fGenParILC = ParToILC(fGenPar);
         /*
-        cout << "ObsTrk::ObsTrk: fGenPar";
-        for (Int_t i = 0; i < 5; i++)cout << fGenPar(i) << ", ";
-        cout << endl;
-        */
-        fObsPar = GenToObsPar(fGenPar, fGC);
-        fObsParMm = ParToMm(fObsPar);
-        fObsParACTS = ParToACTS(fObsPar);
-        fObsParILC = ParToILC(fObsPar);
-        fObsX = ParToX(fObsPar);
-        fObsP = ParToP(fObsPar);
-        fObsQ = ParToQ(fObsPar);
-        fCovMm = CovToMm(fCov);
-        fCovACTS = CovToACTS(fObsPar, fCov);
-        fCovILC = CovToILC(fCov);
-}
-//
-// x[3] track origin, p[3] track momentum at origin, Q charge, B magnetic field in Tesla
-ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC)
-{
-        SetBfield(B);
-        fGC = GC;
-        fGenX.SetXYZ(x[0],x[1],x[2]);
-        fGenP.SetXYZ(p[0],p[1],p[2]);
-        fGenQ = Q;
-        fB = B;
-        fGenPar.ResizeTo(5);
-        fGenParMm.ResizeTo(5);
-        fGenParACTS.ResizeTo(6);
-        fGenParILC.ResizeTo(5);
-        fObsPar.ResizeTo(5);
-        fObsParMm.ResizeTo(5);
-        fObsParACTS.ResizeTo(6);
-        fObsParILC.ResizeTo(5);
-        fCov.ResizeTo(5, 5);
-        fCovMm.ResizeTo(5, 5);
-        fCovACTS.ResizeTo(6, 6);
-        fCovILC.ResizeTo(5, 5);
-        fGenPar = XPtoPar(fGenX, fGenP, Q);
-        fGenParACTS = ParToACTS(fGenPar);
-        fGenParILC = ParToILC(fGenPar);
-        /*
-        cout << "ObsTrk::ObsTrk: fGenPar";
-        for (Int_t i = 0; i < 5; i++)cout << fGenPar(i) << ", ";
-        cout << endl;
+        std::cout << "ObsTrk::ObsTrk: fGenPar";
+        for (Int_t i = 0; i < 5; i++)std::cout << fGenPar(i) << ", ";
+        std::cout << std::endl;
         */
         fObsPar = GenToObsPar(fGenPar, fGC);
@@ -88 +41 @@
         fObsP = ParToP(fObsPar);
         fObsQ = ParToQ(fObsPar);
-        fCovACTS = CovToACTS(fObsPar, fCov);
+        fCovACTS = CovToACTS(fCov);
         fCovILC = CovToILC(fCov);
 }
@@ -98 +51 @@
         fGenP.Clear();
         fGenPar.Clear();
-        fGenParMm.Clear();
         fGenParACTS.Clear();
-        fGenParILC.Clear();
         fObsX.Clear();
         fObsP.Clear();
         fObsPar.Clear();
-        fObsParMm.Clear();
         fObsParACTS.Clear();
-        fObsParILC.Clear();
         fCov.Clear();
-        fCovMm.Clear();
         fCovACTS.Clear();
-        fCovILC.Clear();
+}
+TVectorD ObsTrk::XPtoPar(TVector3 x, TVector3 p, Double_t Q)
+{
+        //
+        TVectorD Par(5);
+        // Transverse parameters
+        Double_t a = -Q*fB*0.2998;                      // Units are Tesla, GeV and meters
+        Double_t pt = p.Pt();
+        Double_t C = a / (2 * pt);                      // Half curvature
+        //std::cout << "ObsTrk::XPtoPar: fB = " << fB << ", a = " << a << ", pt = " << pt << ", C = " << C << std::endl;
+        Double_t r2 = x.Perp2();
+        Double_t cross = x(0)*p(1) - x(1)*p(0);
+        Double_t T = TMath::Sqrt(pt*pt - 2 * a*cross + a*a*r2);
+        Double_t phi0 = TMath::ATan2((p(1) - a*x(0)) / T, (p(0) + a*x(1)) / T); // Phi0
+        Double_t D;                                                     // Impact parameter D
+        if (pt < 10.0) D = (T - pt) / a;
+        else D = (-2 * cross + a*r2) / (T + pt);
+        //
+        Par(0) = D;             // Store D
+        Par(1) = phi0;  // Store phi0
+        Par(2) = C;             // Store C
+        //Longitudinal parameters
+        Double_t B = C*TMath::Sqrt(TMath::Max(r2 - D*D,0.0) / (1 + 2 * C*D));
+        Double_t st = TMath::ASin(B) / C;
+        Double_t ct = p(2) / pt;
+        Double_t z0 = x(2) - ct*st;
+        //
+        Par(3) = z0;            // Store z0
+        Par(4) = ct;            // Store cot(theta)
+        //
+        return Par;
+}
+//
+TVector3 ObsTrk::ParToX(TVectorD Par)
+{
+        Double_t D    = Par(0);
+        Double_t phi0 = Par(1);
+        Double_t z0   = Par(3);
+        //
+        TVector3 Xval;
+        Xval(0) = -D*TMath::Sin(phi0); 
+        Xval(1) =  D*TMath::Cos(phi0);  
+        Xval(2) =  z0;
+        //
+        return Xval;
+}
+//
+TVector3 ObsTrk::ParToP(TVectorD Par)
+{
+        Double_t C    = Par(2);
+        Double_t phi0 = Par(1);
+        Double_t ct   = Par(4);
+        //
+        TVector3 Pval;
+        Double_t pt = fB*0.2998 / TMath::Abs(2 * C);
+        Pval(0) = pt*TMath::Cos(phi0);
+        Pval(1) = pt*TMath::Sin(phi0);
+        Pval(2) = pt*ct;
+        //
+        return Pval;
+}
+//
+
+Double_t ObsTrk::ParToQ(TVectorD Par)
+{
+        return TMath::Sign(1.0, -Par(2));
 }
 //
@@ -122 +135 @@
         // Check ranges
         Double_t minPt = GC->GetMinPt ();
-        //if (pt < minPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << endl;
+        if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;
         Double_t maxPt = GC->GetMaxPt();
-        //if (pt > maxPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << endl;
+        if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;
         Double_t minAn = GC->GetMinAng();
-        //if (angd < minAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd 
-        //      << " is below grid range of " << minAn << endl;
+        if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd 
+                << " is below grid range of " << minAn << std::endl;
         Double_t maxAn = GC->GetMaxAng();
-        //if (angd > maxAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd
-        //      << " is above grid range of " << maxAn << endl;
+        if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
+                << " is above grid range of " << maxAn << std::endl;
         //
         TMatrixDSym Cov = GC->GetCov(pt, angd);
@@ -157 +170 @@
         return oPar;
 }
+// Parameter conversion to ACTS format
+TVectorD ObsTrk::ParToACTS(TVectorD Par)
+{
+        TVectorD pACTS(6);      // Return vector
+        //
+        Double_t b = -0.29988*fB / 2.;
+        pACTS(0) = 1000*Par(0);         // D from m to mm
+        pACTS(1) = 1000 * Par(3);       // z0 from m to mm
+        pACTS(2) = Par(1);                      // Phi0 is unchanged
+        pACTS(3) = TMath::ATan(1.0 / Par(4)) + TMath::PiOver2();                // Theta in [0, pi] range
+        pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4)));         // q/p in GeV
+        pACTS(5) = 0.0;                         // Time: currently undefined
+        //
+        return pACTS;
+}
+// Covariance conversion to ACTS format
+TMatrixDSym ObsTrk::CovToACTS(TMatrixDSym Cov)
+{
+        TMatrixDSym cACTS(6); cACTS.Zero();
+        Double_t b = -0.29988*fB / 2.;
+        //
+        // Fill derivative matrix
+        TMatrixD A(5, 5);       A.Zero();
+        Double_t ct = fGenPar(4);       // cot(theta)
+        Double_t C = fGenPar(2);                // half curvature
+        A(0, 0) = 1000.;                // D-D  conversion to mm
+        A(1, 2) = 1.0;          // phi0-phi0
+        A(2, 4) = 1.0/(TMath::Sqrt(1.0 + ct*ct) * b);   // q/p-C
+        A(3, 1) = 1000.;                // z0-z0 conversion to mm
+        A(4, 3) = -1.0 / (1.0 + ct*ct); // theta - cot(theta)
+        A(4, 4) = -C*ct / (b*pow(1.0 + ct*ct,3.0/2.0)); // q/p-cot(theta)
+        //
+        TMatrixDSym Cv = Cov;
+        TMatrixD At(5, 5);
+        At.Transpose(A);
+        Cv.Similarity(At);
+        TMatrixDSub(cACTS, 0, 4, 0, 4) = Cv;
+        cACTS(5, 5) = 0.1;      // Currently undefined: set to arbitrary value to avoid crashes
+        //
+        return cACTS;
+}
+
+// Parameter conversion to ILC format
+TVectorD ObsTrk::ParToILC(TVectorD Par)
+{
+        TVectorD pILC(5);       // Return vector
+        //
+        pILC(0) = Par(0)*1.0e3;                 // d0 in mm
+        pILC(1) = Par(1);                               // phi0 is unchanged
+        pILC(2) = -2 * Par(2)*1.0e-3;   // w in mm^-1
+        pILC(3) = Par(3)*1.0e3;                 // z0 in mm
+        pILC(4) = Par(4);                               // tan(lambda) = cot(theta)
+        //
+        return pILC;
+}
+// Covariance conversion to ILC format
+TMatrixDSym ObsTrk::CovToILC(TMatrixDSym Cov)
+{
+        TMatrixDSym cILC(5); cILC.Zero();
+        //
+        // Fill derivative matrix
+        TMatrixD A(5, 5);       A.Zero();
+        //
+        A(0, 0) = 1.0e3;                // D-d0 in mm
+        A(1, 1) = 1.0;          // phi0-phi0
+        A(2, 2) = -2.0e-3;      // w-C
+        A(3, 3) = 1.0e3;                // z0-z0 conversion to mm
+        A(4, 4) = 1.0;          // tan(lambda) - cot(theta)
+        //
+        TMatrixDSym Cv = Cov;
+        TMatrixD At(5, 5);
+        At.Transpose(A);
+        Cv.Similarity(At);
+        cILC = Cv;
+        //
+        return cILC;
+}
+
+
+
+