Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • external/TrackCovariance/ObsTrk.cc

    ra0f5d71 rff9fb2d9  
     1#include <iostream>
     2
    13#include <TMath.h>
    24#include <TVector3.h>
    … …  
    57#include <TDecompChol.h>
    68#include <TRandom.h>
    7 #include <iostream>
     9
    810#include "SolGridCov.h"
    911#include "ObsTrk.h"
    10 //
    11 // Constructors
     12
     13using namespace std;
     14
    1215// x(3) track origin, p(3) track momentum at origin, Q charge, B magnetic field in Tesla
    1316ObsTrk::ObsTrk(TVector3 x, TVector3 p, Double_t Q, Double_t B, SolGridCov *GC)
    1417{
    15         fGC = GC;
    16         fGenX = x;
    17         fGenP = p;
    18         fGenQ = Q;
    19         fB = B;
    20         fGenPar.ResizeTo(5);
    21         fGenParACTS.ResizeTo(6);
    22         fGenParILC.ResizeTo(5);
    23         fObsPar.ResizeTo(5);
    24         fObsParACTS.ResizeTo(6);
    25         fObsParILC.ResizeTo(5);
    26         fCov.ResizeTo(5, 5);
    27         fCovACTS.ResizeTo(6, 6);
    28         fCovILC.ResizeTo(5, 5);
    29         fGenPar = XPtoPar(x,p,Q);
    30         fGenParACTS = ParToACTS(fGenPar);
    31         fGenParILC = ParToILC(fGenPar);
    32         /*
    33         std::cout << "ObsTrk::ObsTrk: fGenPar";
    34         for (Int_t i = 0; i < 5; i++)std::cout << fGenPar(i) << ", ";
    35         std::cout << std::endl;
    36         */
    37         fObsPar = GenToObsPar(fGenPar, fGC);
    38         fObsParACTS = ParToACTS(fObsPar);
    39         fObsParILC = ParToILC(fObsPar);
    40         fObsX = ParToX(fObsPar);
    41         fObsP = ParToP(fObsPar);
    42         fObsQ = ParToQ(fObsPar);
    43         fCovACTS = CovToACTS(fCov);
    44         fCovILC = CovToILC(fCov);
     18  fGC = GC;
     19  fGenX = x;
     20  fGenP = p;
     21  fGenQ = Q;
     22  fB = B;
     23  fGenPar.ResizeTo(5);
     24  fObsPar.ResizeTo(5);
     25  fCov.ResizeTo(5, 5);
     26  fGenPar = XPtoPar(x, p, Q);
     27  fObsPar = GenToObsPar(fGenPar, fGC);
     28  fObsX = ParToX(fObsPar);
     29  fObsP = ParToP(fObsPar);
     30  fObsQ = ParToQ(fObsPar);
    4531}
    46 //
    47 // Destructor
     32
    4833ObsTrk::~ObsTrk()
    4934{
    50         fGenX.Clear();
    51         fGenP.Clear();
    52         fGenPar.Clear();
    53         fGenParACTS.Clear();
    54         fObsX.Clear();
    55         fObsP.Clear();
    56         fObsPar.Clear();
    57         fObsParACTS.Clear();
    58         fCov.Clear();
    59         fCovACTS.Clear();
    6035}
     36
    6137TVectorD ObsTrk::XPtoPar(TVector3 x, TVector3 p, Double_t Q)
    6238{
    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;
     39  TVectorD Par(5);
     40  // Transverse parameters
     41  Double_t a = -Q * fB * 0.2998; // Units are Tesla, GeV and m
     42  Double_t pt = p.Pt();
     43  Double_t C = a / (2 * pt); // Half curvature
     44
     45  Double_t r2 = x.Perp2();
     46  Double_t cross = x(0) * p(1) - x(1) * p(0);
     47  Double_t T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2);
     48  Double_t phi0 = TMath::ATan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); // Phi0
     49  Double_t D; // Impact parameter D
     50  if (pt < 10.0) D = (T - pt) / a;
     51  else D = (-2 * cross + a * r2) / (T + pt);
     52
     53  Par(0) = D; // Store D
     54  Par(1) = phi0; // Store phi0
     55  Par(2) = C; // Store C
     56  // Longitudinal parameters
     57  Double_t B = C * TMath::Sqrt(TMath::Max(r2 - D * D,0.0) / (1 + 2 * C * D));
     58  Double_t st = TMath::ASin(B) / C;
     59  Double_t ct = p(2) / pt;
     60  Double_t z0 = x(2) - ct * st;
     61
     62  Par(3) = z0; // Store z0
     63  Par(4) = ct; // Store cot(theta)
     64
     65  return Par;
    9166}
    92 //
     67
    9368TVector3 ObsTrk::ParToX(TVectorD Par)
    9469{
    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;
     70  Double_t D    = Par(0);
     71  Double_t phi0 = Par(1);
     72  Double_t z0   = Par(3);
     73
     74  TVector3 Xval;
     75  Xval(0) = -D * TMath::Sin(phi0);
     76  Xval(1) =  D * TMath::Cos(phi0);
     77  Xval(2) =  z0;
     78
     79  return Xval;
    10580}
    106 //
     81
    10782TVector3 ObsTrk::ParToP(TVectorD Par)
    10883{
    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;
     84  Double_t C    = Par(2);
     85  Double_t phi0 = Par(1);
     86  Double_t ct   = Par(4);
     87  //
     88  TVector3 Pval;
     89  Double_t pt = fB * 0.2998 / TMath::Abs(2 * C);
     90  Pval(0) = pt * TMath::Cos(phi0);
     91  Pval(1) = pt * TMath::Sin(phi0);
     92  Pval(2) = pt * ct;
     93  //
     94  return Pval;
    12095}
    121 //
    12296
    12397Double_t ObsTrk::ParToQ(TVectorD Par)
    12498{
    125         return TMath::Sign(1.0, -Par(2));
     99  return TMath::Sign(1.0, -Par(2));
    126100}
    127 //
     101
    128102TVectorD ObsTrk::GenToObsPar(TVectorD gPar, SolGridCov *GC)
    129103{
    130         TVector3 p = ParToP(gPar);
    131         Double_t pt = p.Pt();
    132         Double_t tanTh = 1.0 / TMath::Abs(gPar(4));
    133         Double_t angd = TMath::ATan(tanTh)*180. / TMath::Pi();
    134         //
    135         // Check ranges
    136         Double_t minPt = GC->GetMinPt ();
    137         if (pt < minPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << std::endl;
    138         Double_t maxPt = GC->GetMaxPt();
    139         if (pt > maxPt) std::cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << std::endl;
    140         Double_t minAn = GC->GetMinAng();
    141         if (angd < minAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
    142                 << " is below grid range of " << minAn << std::endl;
    143         Double_t maxAn = GC->GetMaxAng();
    144         if (angd > maxAn) std::cout << "Warning ObsTrk::GenToObsPar: angle " << angd
    145                 << " is above grid range of " << maxAn << std::endl;
    146         //
    147         TMatrixDSym Cov = GC->GetCov(pt, angd);
    148         fCov = Cov;
    149         //
    150         // Now do Choleski decomposition and random number extraction, with appropriate stabilization
    151         //
    152         TMatrixDSym CvN = Cov;
    153         TMatrixDSym DCv(5); DCv.Zero();
    154         TMatrixDSym DCvInv(5); DCvInv.Zero();
    155         for (Int_t id = 0; id < 5; id++)
    156         {
    157                 Double_t dVal = TMath::Sqrt(Cov(id, id));
    158                 DCv   (id, id) = dVal;
    159                 DCvInv(id, id) = 1.0 / dVal;
    160         }
    161         CvN.Similarity(DCvInv);                 // Normalize diagonal to 1
    162         TDecompChol Chl(CvN);
    163         Bool_t OK = Chl.Decompose();            // Choleski decomposition of normalized matrix
    164         TMatrixD U = Chl.GetU();                        // Get Upper triangular matrix
    165         TMatrixD Ut(TMatrixD::kTransposed, U); // Transposed of U (lower triangular)
    166         TVectorD r(5);
    167         for (Int_t i = 0; i < 5; i++)r(i) = gRandom->Gaus(0.0, 1.0);            // Array of normal random numbers
    168         TVectorD oPar = gPar + DCv*(Ut*r);      // Observed parameter vector
    169         //
    170         return oPar;
     104  TVector3 p = ParToP(gPar);
     105  Double_t pt = p.Pt();
     106  Double_t tanTh = 1.0 / TMath::Abs(gPar(4));
     107  Double_t angd = TMath::ATan(tanTh) * 180. / TMath::Pi();
     108  // Check ranges
     109  Double_t minPt = GC->GetMinPt ();
     110  if (pt < minPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is below grid range of " << minPt << endl;
     111  Double_t maxPt = GC->GetMaxPt();
     112  if (pt > maxPt) cout << "Warning ObsTrk::GenToObsPar: pt " << pt << " is above grid range of " << maxPt << endl;
     113  Double_t minAn = GC->GetMinAng();
     114  if (angd < minAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd
     115    << " is below grid range of " << minAn << endl;
     116  Double_t maxAn = GC->GetMaxAng();
     117  if (angd > maxAn) cout << "Warning ObsTrk::GenToObsPar: angle " << angd
     118    << " is above grid range of " << maxAn << endl;
     119  TMatrixDSym Cov = GC->GetCov(pt, angd);
     120  fCov = Cov;
     121  // Now do Choleski decomposition and random number extraction, with appropriate stabilization
     122  TMatrixDSym CvN = Cov;
     123  TMatrixDSym DCv(5); DCv.Zero();
     124  TMatrixDSym DCvInv(5); DCvInv.Zero();
     125  for (Int_t id = 0; id < 5; id++)
     126  {
     127    Double_t dVal = TMath::Sqrt(Cov(id, id));
     128    DCv   (id, id) = dVal;
     129    DCvInv(id, id) = 1.0 / dVal;
     130  }
     131  CvN.Similarity(DCvInv); // Normalize diagonal to 1
     132  TDecompChol Chl(CvN);
     133  Bool_t OK = Chl.Decompose(); // Choleski decomposition of normalized matrix
     134  TMatrixD U = Chl.GetU(); // Get Upper triangular matrix
     135  TMatrixD Ut(TMatrixD::kTransposed, U); // Transposed of U (lower triangular)
     136  TVectorD r(5);
     137  for (Int_t i = 0; i < 5; i++) r(i) = gRandom->Gaus(0.0, 1.0); // Array of normal random numbers
     138  TVectorD oPar = gPar + DCv * (Ut * r); // Observed parameter vector
     139
     140  return oPar;
    171141}
    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.