@@ -527,8 +527,9 @@ GPUd() void TrackParametrizationWithError<value_T>::set(const dim3_t& xyz, const
527527 math_utils::detail::rotateZ<value_t >(ver, -alp);
528528 math_utils::detail::rotateZ<value_t >(mom, -alp);
529529 //
530- value_t pt = gpu::CAMath::Sqrt (mom[0 ] * mom[0 ] + mom[1 ] * mom[1 ]);
531- value_t ptI = 1 .f / pt;
530+ const value_t pt2 = mom[0 ] * mom[0 ] + mom[1 ] * mom[1 ];
531+ const value_t pt = gpu::CAMath::Sqrt (pt2);
532+ const value_t ptI = 1 .f / pt;
532533 this ->setX (ver[0 ]);
533534 this ->setAlpha (alp);
534535 this ->setY (ver[1 ]);
@@ -545,89 +546,53 @@ GPUd() void TrackParametrizationWithError<value_T>::set(const dim3_t& xyz, const
545546 this ->setSnp (-1 .f + kSafe ); // Protection
546547 }
547548 //
548- // Covariance matrix (formulas to be simplified)
549- value_t r = mom[0 ] * ptI; // cos(phi)
550- value_t cv34 = gpu::CAMath::Sqrt (cv[3 ] * cv[3 ] + cv[4 ] * cv[4 ]);
551- //
552- int special = 0 ;
553- value_t sgcheck = r * sn + this ->getSnp () * cs;
554- if (gpu::CAMath::Abs (sgcheck) > 1 - kSafe ) { // special case: lab phi is +-pi/2
555- special = 1 ;
556- sgcheck = sgcheck < 0 ? -1 .f : 1 .f ;
557- } else if (gpu::CAMath::Abs (sgcheck) < kSafe ) {
558- sgcheck = cs < 0 ? -1 .0f : 1 .0f ;
559- special = 2 ; // special case: lab phi is 0
549+ // Covariance matrix from the fixed-alpha Jacobian
550+ // d(Y,Z,snp,tgl,q/pt) / d(X,Y,Z,Px,Py,Pz).
551+ const value_t pt3I = ptI / pt2;
552+ const value_t qeff = charge ? static_cast <value_t >(charge) : 1 .f ;
553+
554+ value_t cLab[6 ][6 ] = {};
555+ int idx = 0 ;
556+ for (int i = 0 ; i < 6 ; ++i) {
557+ for (int j = 0 ; j <= i; ++j) {
558+ cLab[i][j] = cLab[j][i] = cv[idx++];
559+ }
560560 }
561- //
562- mC [kSigY2 ] = cv[0 ] + cv[2 ];
563- mC [kSigZY ] = (-cv[3 ] * sn) < 0 ? -cv34 : cv34;
564- mC [kSigZ2 ] = cv[5 ];
565- //
566- value_t ptI2 = ptI * ptI;
567- value_t tgl2 = this ->getTgl () * this ->getTgl ();
568- if (special == 1 ) {
569- mC [kSigSnpY ] = cv[6 ] * ptI;
570- mC [kSigSnpZ ] = -sgcheck * cv[8 ] * r * ptI;
571- mC [kSigSnp2 ] = gpu::CAMath::Abs (cv[9 ] * r * r * ptI2);
572- mC [kSigTglY ] = (cv[10 ] * this ->getTgl () - sgcheck * cv[15 ]) * ptI / r;
573- mC [kSigTglZ ] = (cv[17 ] - sgcheck * cv[12 ] * this ->getTgl ()) * ptI;
574- mC [kSigTglSnp ] = (-sgcheck * cv[18 ] + cv[13 ] * this ->getTgl ()) * r * ptI2;
575- mC [kSigTgl2 ] = gpu::CAMath::Abs (cv[20 ] - 2 * sgcheck * cv[19 ] * mC [4 ] + cv[14 ] * tgl2) * ptI2;
576- mC [kSigQ2PtY ] = cv[10 ] * ptI2 / r * charge;
577- mC [kSigQ2PtZ ] = -sgcheck * cv[12 ] * ptI2 * charge;
578- mC [kSigQ2PtSnp ] = cv[13 ] * r * ptI * ptI2 * charge;
579- mC [kSigQ2PtTgl ] = (-sgcheck * cv[19 ] + cv[14 ] * this ->getTgl ()) * r * ptI2 * ptI;
580- mC [kSigQ2Pt2 ] = gpu::CAMath::Abs (cv[14 ] * ptI2 * ptI2);
581- } else if (special == 2 ) {
582- mC [kSigSnpY ] = -cv[10 ] * ptI * cs / sn;
583- mC [kSigSnpZ ] = cv[12 ] * cs * ptI;
584- mC [kSigSnp2 ] = gpu::CAMath::Abs (cv[14 ] * cs * cs * ptI2);
585- mC [kSigTglY ] = (sgcheck * cv[6 ] * this ->getTgl () - cv[15 ]) * ptI / sn;
586- mC [kSigTglZ ] = (cv[17 ] - sgcheck * cv[8 ] * this ->getTgl ()) * ptI;
587- mC [kSigTglSnp ] = (cv[19 ] - sgcheck * cv[13 ] * this ->getTgl ()) * cs * ptI2;
588- mC [kSigTgl2 ] = gpu::CAMath::Abs (cv[20 ] - 2 * sgcheck * cv[18 ] * this ->getTgl () + cv[9 ] * tgl2) * ptI2;
589- mC [kSigQ2PtY ] = sgcheck * cv[6 ] * ptI2 / sn * charge;
590- mC [kSigQ2PtZ ] = -sgcheck * cv[8 ] * ptI2 * charge;
591- mC [kSigQ2PtSnp ] = -sgcheck * cv[13 ] * cs * ptI * ptI2 * charge;
592- mC [kSigQ2PtTgl ] = (-sgcheck * cv[18 ] + cv[9 ] * this ->getTgl ()) * ptI2 * ptI * charge;
593- mC [kSigQ2Pt2 ] = gpu::CAMath::Abs (cv[9 ] * ptI2 * ptI2);
594- } else {
595- double m00 = -sn; // m10=cs;
596- double m23 = -pt * (sn + this ->getSnp () * cs / r), m43 = -pt * pt * (r * cs - this ->getSnp () * sn);
597- double m24 = pt * (cs - this ->getSnp () * sn / r), m44 = -pt * pt * (r * sn + this ->getSnp () * cs);
598- double m35 = pt, m45 = -pt * pt * this ->getTgl ();
599- //
600- if (charge) { // RS: this is a hack, proper treatment to be implemented
601- m43 *= charge;
602- m44 *= charge;
603- m45 *= charge;
561+
562+ value_t jac[5 ][6 ] = {};
563+ jac[kY ][0 ] = -sn;
564+ jac[kY ][1 ] = cs;
565+ jac[kZ ][2 ] = 1 .;
566+
567+ const value_t u = mom[0 ];
568+ const value_t v = mom[1 ];
569+ const value_t w = mom[2 ];
570+ const value_t dSnpDu = -u * v * pt3I;
571+ const value_t dSnpDv = u * u * pt3I;
572+ const value_t dTglDu = -w * u * pt3I;
573+ const value_t dTglDv = -w * v * pt3I;
574+ const value_t dTglDw = ptI;
575+ const value_t dQ2PtDu = -qeff * u * pt3I;
576+ const value_t dQ2PtDv = -qeff * v * pt3I;
577+
578+ jac[kSnp ][3 ] = dSnpDu * cs - dSnpDv * sn;
579+ jac[kSnp ][4 ] = dSnpDu * sn + dSnpDv * cs;
580+ jac[kTgl ][3 ] = dTglDu * cs - dTglDv * sn;
581+ jac[kTgl ][4 ] = dTglDu * sn + dTglDv * cs;
582+ jac[kTgl ][5 ] = dTglDw;
583+ jac[kQ2Pt ][3 ] = dQ2PtDu * cs - dQ2PtDv * sn;
584+ jac[kQ2Pt ][4 ] = dQ2PtDu * sn + dQ2PtDv * cs;
585+
586+ for (int i = 0 ; i < kNParams ; ++i) {
587+ for (int j = 0 ; j <= i; ++j) {
588+ value_t cij = 0 .;
589+ for (int k = 0 ; k < 6 ; ++k) {
590+ for (int l = 0 ; l < 6 ; ++l) {
591+ cij += jac[i][k] * cLab[k][l] * jac[j][l];
592+ }
593+ }
594+ mC [CovarMap[i][j]] = cij;
604595 }
605- //
606- double a1 = cv[13 ] - cv[9 ] * (m23 * m44 + m43 * m24) / m23 / m43;
607- double a2 = m23 * m24 - m23 * (m23 * m44 + m43 * m24) / m43;
608- double a3 = m43 * m44 - m43 * (m23 * m44 + m43 * m24) / m23;
609- double a4 = cv[14 ] + 2 . * cv[9 ];
610- double a5 = m24 * m24 - 2 . * m24 * m44 * m23 / m43;
611- double a6 = m44 * m44 - 2 . * m24 * m44 * m43 / m23;
612- //
613- mC [kSigSnpY ] = (cv[10 ] * m43 - cv[6 ] * m44) / (m24 * m43 - m23 * m44) / m00;
614- mC [kSigQ2PtY ] = (cv[6 ] / m00 - mC [kSigSnpY ] * m23) / m43;
615- mC [kSigTglY ] = (cv[15 ] / m00 - mC [kSigQ2PtY ] * m45) / m35;
616- mC [kSigSnpZ ] = (cv[12 ] * m43 - cv[8 ] * m44) / (m24 * m43 - m23 * m44);
617- mC [kSigQ2PtZ ] = (cv[8 ] - mC [kSigSnpZ ] * m23) / m43;
618- mC [kSigTglZ ] = cv[17 ] / m35 - mC [kSigQ2PtZ ] * m45 / m35;
619- mC [kSigSnp2 ] = gpu::CAMath::Abs ((a4 * a3 - a6 * a1) / (a5 * a3 - a6 * a2));
620- mC [kSigQ2Pt2 ] = gpu::CAMath::Abs ((a1 - a2 * mC [kSigSnp2 ]) / a3);
621- mC [kSigQ2PtSnp ] = (cv[9 ] - mC [kSigSnp2 ] * m23 * m23 - mC [kSigQ2Pt2 ] * m43 * m43) / m23 / m43;
622- double b1 = cv[18 ] - mC [kSigQ2PtSnp ] * m23 * m45 - mC [kSigQ2Pt2 ] * m43 * m45;
623- double b2 = m23 * m35;
624- double b3 = m43 * m35;
625- double b4 = cv[19 ] - mC [kSigQ2PtSnp ] * m24 * m45 - mC [kSigQ2Pt2 ] * m44 * m45;
626- double b5 = m24 * m35;
627- double b6 = m44 * m35;
628- mC [kSigTglSnp ] = (b4 - b6 * b1 / b3) / (b5 - b6 * b2 / b3);
629- mC [kSigQ2PtTgl ] = b1 / b3 - b2 * mC [kSigTglSnp ] / b3;
630- mC [kSigTgl2 ] = gpu::CAMath::Abs ((cv[20 ] - mC [kSigQ2Pt2 ] * (m45 * m45) - mC [kSigQ2PtTgl ] * 2 .f * m35 * m45) / (m35 * m35));
631596 }
632597 checkCovariance ();
633598}
@@ -1668,42 +1633,52 @@ GPUd() bool TrackParametrizationWithError<value_T>::getCovXYZPxPyPzGlo(std::arra
16681633 return false ;
16691634 }
16701635
1671- auto pt = this ->getPt ();
1672- value_t sn, cs;
1636+ const value_t pt = this ->getPt ();
1637+ const value_t q2pt = this ->getQ2Pt ();
1638+ value_t sn = 0 .f , cs = 0 .f ;
16731639 o2::math_utils::detail::sincos (this ->getAlpha (), sn, cs);
1674- auto r = gpu::CAMath::Sqrt ((1 . - this ->getSnp ()) * (1 . + this ->getSnp ()));
1675- auto m00 = -sn, m10 = cs;
1676- auto m23 = -pt * (sn + this ->getSnp () * cs / r), m43 = -pt * pt * (r * cs - this ->getSnp () * sn);
1677- auto m24 = pt * (cs - this ->getSnp () * sn / r), m44 = -pt * pt * (r * sn + this ->getSnp () * cs);
1678- auto m35 = pt, m45 = -pt * pt * this ->getTgl ();
1679-
1680- if (this ->getSign () < 0 ) {
1681- m43 = -m43;
1682- m44 = -m44;
1683- m45 = -m45;
1684- }
1685-
1686- cv[0 ] = mC [0 ] * m00 * m00;
1687- cv[1 ] = mC [0 ] * m00 * m10;
1688- cv[2 ] = mC [0 ] * m10 * m10;
1689- cv[3 ] = mC [1 ] * m00;
1690- cv[4 ] = mC [1 ] * m10;
1691- cv[5 ] = mC [2 ];
1692- cv[6 ] = m00 * (mC [3 ] * m23 + mC [10 ] * m43);
1693- cv[7 ] = m10 * (mC [3 ] * m23 + mC [10 ] * m43);
1694- cv[8 ] = mC [4 ] * m23 + mC [11 ] * m43;
1695- cv[9 ] = m23 * (mC [5 ] * m23 + mC [12 ] * m43) + m43 * (mC [12 ] * m23 + mC [14 ] * m43);
1696- cv[10 ] = m00 * (mC [3 ] * m24 + mC [10 ] * m44);
1697- cv[11 ] = m10 * (mC [3 ] * m24 + mC [10 ] * m44);
1698- cv[12 ] = mC [4 ] * m24 + mC [11 ] * m44;
1699- cv[13 ] = m23 * (mC [5 ] * m24 + mC [12 ] * m44) + m43 * (mC [12 ] * m24 + mC [14 ] * m44);
1700- cv[14 ] = m24 * (mC [5 ] * m24 + mC [12 ] * m44) + m44 * (mC [12 ] * m24 + mC [14 ] * m44);
1701- cv[15 ] = m00 * (mC [6 ] * m35 + mC [10 ] * m45);
1702- cv[16 ] = m10 * (mC [6 ] * m35 + mC [10 ] * m45);
1703- cv[17 ] = mC [7 ] * m35 + mC [11 ] * m45;
1704- cv[18 ] = m23 * (mC [8 ] * m35 + mC [12 ] * m45) + m43 * (mC [13 ] * m35 + mC [14 ] * m45);
1705- cv[19 ] = m24 * (mC [8 ] * m35 + mC [12 ] * m45) + m44 * (mC [13 ] * m35 + mC [14 ] * m45);
1706- cv[20 ] = m35 * (mC [9 ] * m35 + mC [13 ] * m45) + m45 * (mC [13 ] * m35 + mC [14 ] * m45);
1640+ const value_t snp = this ->getSnp ();
1641+ const value_t csp = gpu::CAMath::Sqrt ((1 .f - snp) * (1 .f + snp));
1642+ const value_t pXLoc = pt * csp;
1643+ const value_t pYLoc = pt * snp;
1644+ const value_t pZ = pt * this ->getTgl ();
1645+ const value_t pX = cs * pXLoc - sn * pYLoc;
1646+ const value_t pY = sn * pXLoc + cs * pYLoc;
1647+
1648+ value_t cTr[5 ][5 ] = {};
1649+ for (int i = 0 ; i < kNParams ; ++i) {
1650+ for (int j = 0 ; j <= i; ++j) {
1651+ cTr[i][j] = cTr[j][i] = mC [CovarMap[i][j]];
1652+ }
1653+ }
1654+
1655+ double jac[6 ][5 ] = {};
1656+ jac[0 ][kY ] = -sn;
1657+ jac[1 ][kY ] = cs;
1658+ jac[2 ][kZ ] = 1 .f ;
1659+
1660+ const value_t dPxDSnp = -pt * (cs * snp / csp + sn);
1661+ const value_t dPyDSnp = pt * (cs - sn * snp / csp);
1662+ jac[3 ][kSnp ] = dPxDSnp;
1663+ jac[4 ][kSnp ] = dPyDSnp;
1664+ jac[5 ][kTgl ] = pt;
1665+
1666+ jac[3 ][kQ2Pt ] = -pX / q2pt;
1667+ jac[4 ][kQ2Pt ] = -pY / q2pt;
1668+ jac[5 ][kQ2Pt ] = -pZ / q2pt;
1669+
1670+ int idx = 0 ;
1671+ for (int i = 0 ; i < 6 ; ++i) {
1672+ for (int j = 0 ; j <= i; ++j) {
1673+ double cij = 0 .f ;
1674+ for (int k = 0 ; k < kNParams ; ++k) {
1675+ for (int l = 0 ; l < kNParams ; ++l) {
1676+ cij += jac[i][k] * cTr[k][l] * jac[j][l];
1677+ }
1678+ }
1679+ cv[idx++] = cij;
1680+ }
1681+ }
17071682
17081683 return true ;
17091684}
0 commit comments