Fork me on GitHub

Changeset a617744 in git for external/TrackCovariance


Ignore:
Timestamp:
Mar 1, 2021, 4:01:37 PM (4 years ago)
Author:
michele <michele.selvaggi@…>
Branches:
master
Children:
40e890c
Parents:
3f8466a
Message:

adding latest TrackCovariance libraries (F. Bedeschi)

Location:
external/TrackCovariance
Files:
2 added
2 edited

Legend:

Unmodified
Added
Removed

  • external/TrackCovariance/ObsTrk.cc

    r3f8466a ra617744  
    1010//
    1111// Constructors
     12//
    1213// x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla
    1314ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC)
    1415{
     16        SetBfield(B);
    1517        fGC = GC;
    1618        fGenX = x;
    … …  
    1921        fB = B;
    2022        fGenPar.ResizeTo(5);
     23        fGenParMm.ResizeTo(5);
    2124        fGenParACTS.ResizeTo(6);
    2225        fGenParILC.ResizeTo(5);
    2326        fObsPar.ResizeTo(5);
     27        fObsParMm.ResizeTo(5);
    2428        fObsParACTS.ResizeTo(6);
    2529        fObsParILC.ResizeTo(5);
    2630        fCov.ResizeTo(5, 5);
     31        fCovMm.ResizeTo(5, 5);
    2732        fCovACTS.ResizeTo(6, 6);
    2833        fCovILC.ResizeTo(5, 5);
    2934        fGenPar = XPtoPar(x,p,Q);
     35        fGenParMm = ParToMm(fGenPar);
    3036        fGenParACTS = ParToACTS(fGenPar);
    3137        fGenParILC = ParToILC(fGenPar);
    3238        /*
    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;
    3642        */
    3743        fObsPar = GenToObsPar(fGenPar, fGC);
     44        fObsParMm = ParToMm(fObsPar);
    3845        fObsParACTS = ParToACTS(fObsPar);
    3946        fObsParILC = ParToILC(fObsPar);
    … …  
    4148        fObsP = ParToP(fObsPar);
    4249        fObsQ = ParToQ(fObsPar);
    43         fCovACTS = CovToACTS(fCov);
     50        fCovMm = CovToMm(fCov);
     51        fCovACTS = CovToACTS(fObsPar, fCov);
    4452        fCovILC = CovToILC(fCov);
    4553}
    46 
     54//
    4755// x[3] track origin, p[3] track momentum at origin, Q charge, B magnetic field in Tesla
    4856ObsTrk::ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC)
    4957{
     58        SetBfield(B);
    5059        fGC = GC;
    5160        fGenX.SetXYZ(x[0],x[1],x[2]);
    … …  
    5463        fB = B;
    5564        fGenPar.ResizeTo(5);
     65        fGenParMm.ResizeTo(5);
    5666        fGenParACTS.ResizeTo(6);
    5767        fGenParILC.ResizeTo(5);
    5868        fObsPar.ResizeTo(5);
     69        fObsParMm.ResizeTo(5);
    5970        fObsParACTS.ResizeTo(6);
    6071        fObsParILC.ResizeTo(5);
    6172        fCov.ResizeTo(5, 5);
     73        fCovMm.ResizeTo(5, 5);
    6274        fCovACTS.ResizeTo(6, 6);
    6375        fCovILC.ResizeTo(5, 5);
    … …  
    7688        fObsP = ParToP(fObsPar);
    7789        fObsQ = ParToQ(fObsPar);
    78         fCovACTS = CovToACTS(fCov);
     90        fCovACTS = CovToACTS(fObsPar, fCov);
    7991        fCovILC = CovToILC(fCov);
    8092}
    81 
    82 
    8393//
    8494// Destructor
    … …  
    8898        fGenP.Clear();
    8999        fGenPar.Clear();
     100        fGenParMm.Clear();
    90101        fGenParACTS.Clear();
     102        fGenParILC.Clear();
    91103        fObsX.Clear();
    92104        fObsP.Clear();
    93105        fObsPar.Clear();
     106        fObsParMm.Clear();
    94107        fObsParACTS.Clear();
     108        fObsParILC.Clear();
    95109        fCov.Clear();
     110        fCovMm.Clear();
    96111        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();
    163113}
    164114//
    … …  
    172122        // Check ranges
    173123        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;
    175125        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;
    177127        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;
    180130        Double_t maxAn = GC->GetMaxAng();
    181         //if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
    182         //      << " 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;
    183133        //
    184134        TMatrixDSym Cov = GC->GetCov(pt, angd);
    … …  
    207157        return oPar;
    208158}
    209 // Parameter conversion to ACTS format
    210 TVectorD ObsTrk::ParToACTS(TVectorD Par)
    211 {
    212         TVectorD pACTS(6);      // Return vector
    213         //
    214         Double_t b = -0.29988*fB / 2.;
    215         pACTS(0) = 1000*Par(0);         // D from m to mm
    216         pACTS(1) = 1000 * Par(3);       // z0 from m to mm
    217         pACTS(2) = Par(1);                      // Phi0 is unchanged
    218         pACTS(3) = TMath::ATan2(1.0,Par(4));            // Theta in [0, pi] range
    219         pACTS(4) = Par(2) / (b*TMath::Sqrt(1 + Par(4)*Par(4)));         // q/p in GeV
    220         pACTS(5) = 0.0;                         // Time: currently undefined
    221         //
    222         return pACTS;
    223 }
    224 // Covariance conversion to ACTS format
    225 TMatrixDSym ObsTrk::CovToACTS(TMatrixDSym Cov)
    226 {
    227         TMatrixDSym cACTS(6); cACTS.Zero();
    228         Double_t b = -0.29988*fB / 2.;
    229         //
    230         // Fill derivative matrix
    231         TMatrixD A(5, 5);       A.Zero();
    232         Double_t ct = fGenPar(4);       // cot(theta)
    233         Double_t C = fGenPar(2);                // half curvature
    234         A(0, 0) = 1000.;                // D-D  conversion to mm
    235         A(1, 2) = 1.0;          // phi0-phi0
    236         A(2, 4) = 1.0/(TMath::Sqrt(1.0 + ct*ct) * b);   // q/p-C
    237         A(3, 1) = 1000.;                // z0-z0 conversion to mm
    238         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 crashes
    247         //
    248         return cACTS;
    249 }
    250 
    251 // Parameter conversion to ILC format
    252 TVectorD ObsTrk::ParToILC(TVectorD Par)
    253 {
    254         TVectorD pILC(5);       // Return vector
    255         //
    256         pILC(0) = Par(0)*1.0e3;                 // d0 in mm
    257         pILC(1) = Par(1);                               // phi0 is unchanged
    258         pILC(2) = -2 * Par(2)*1.0e-3;   // w in mm^-1
    259         pILC(3) = Par(3)*1.0e3;                 // z0 in mm
    260         pILC(4) = Par(4);                               // tan(lambda) = cot(theta)
    261         //
    262         return pILC;
    263 }
    264 // Covariance conversion to ILC format
    265 TMatrixDSym ObsTrk::CovToILC(TMatrixDSym Cov)
    266 {
    267         TMatrixDSym cILC(5); cILC.Zero();
    268         //
    269         // Fill derivative matrix
    270         TMatrixD A(5, 5);       A.Zero();
    271         //
    272         A(0, 0) = 1.0e3;                // D-d0 in mm
    273         A(1, 1) = 1.0;          // phi0-phi0
    274         A(2, 2) = -2.0e-3;      // w-C
    275         A(3, 3) = 1.0e3;                // z0-z0 conversion to mm
    276         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  
    66#include <TMatrixDSym.h>
    77#include <TDecompChol.h>
     8#include "TrkUtil.h"
    89#include "SolGridCov.h"
    910//
    … …  
    1314//         INFN - Sezione di Pisa, Italy
    1415//
    15 class ObsTrk{
     16class ObsTrk: public TrkUtil
     17{
    1618        //
    1719        // Class to handle simulation of tracking resolution
    … …  
    1921        // Prefix Gen marks variables before resolution smearing
    2022        //
    21 private:
    22         Double_t fB;                                            // Solenoid magnetic field
    23         SolGridCov *fGC;                                        // Covariance matrix grid
     23private:       
     24        Double_t fB;                                    // Solenoid magnetic field
     25        SolGridCov *fGC;                                // Covariance matrix grid
    2426        Double_t fGenQ;                                 // Generated track charge
    2527        Double_t fObsQ;                                 // Observed  track charge
    2628        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 
    2830        TVector3 fGenP;                                 // Generated track momentum at track origin
    2931        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
    3134        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
    3438        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
    3742        TMatrixDSym fCovACTS;                   // Covariance of track parameters in ACTS format
    3843                                                                        // (D, z0, phi0, theta, q/p, time)
    39         TMatrixDSym fCovILC;                            // Covariance of track parameters in ILC format
     44        TMatrixDSym fCovILC;                    // Covariance of track parameters in ILC format
    4045                                                                        // (d0, phi0, w, z0, tan(lambda))
    4146        //
    42         // Conversion to ACTS parametrization
     47        // Service routines
    4348        //
    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);
    5150        //
    5251public:
    … …  
    5453        // Constructors
    5554        // 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 track
     55        ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC); // Initialize and generate smeared
    5756        ObsTrk(Double_t *x, Double_t *p, Double_t Q, Double_t B, SolGridCov* GC);       // Initialize and generate smeared track
    58   // Destructor
     57        // Destructor
    5958        ~ObsTrk();
    60         //
    61         // Service routines
    62         //
    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);
    6859        //
    6960        // Accessors
    … …  
    7566        TVector3 GetGenP()      { return fGenP; }
    7667        // 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
    7870        // D, z0, phi0, theta, q/p, time
    7971        TVectorD GetGenParACTS()        { return fGenParACTS; }
    … …  
    8577        TVector3 GetObsP()      { return fObsP; }
    8678        // 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
    8881        // D, z0, phi0, theta, q/p, time
    8982        TVectorD GetObsParACTS()        { return fObsParACTS; }
    … …  
    9184        TVectorD GetObsParILC() { return fObsParILC; }
    9285        // Covariances
    93         TMatrixDSym GetCov(){ return fCov; }
     86        TMatrixDSym GetCov()    { return fCov; }        // in meters
     87        TMatrixDSym GetCovMm()  { return fCov; }        // in mm
    9488        TMatrixDSym GetCovACTS(){ return fCovACTS; }
    9589        TMatrixDSym GetCovILC(){ return fCovILC; }
     90        //
    9691};
    9792
Note: See TracChangeset for help on using the changeset viewer.