Skip to content
Snippets Groups Projects

Add JacobdP4dMom version considering charge

Merged Hendrik Jage requested to merge hjage_add_ions into 2024-patches
@@ -632,6 +632,74 @@ namespace Gaudi {
return;
}
// ========================================================================
/** Compute the jacobian for the transformation of a covariance matrix
* with rows representing track parameters TxTyQop and columns in xyz
* into a covariance matrix representing the track parameters in
* PxPyPzE and columns xyz, while taking the particles charge into
* account. Based on equivalent without charge from Sean BRISBANE.
*
* @code
*
* ROOT::Math::SMatrix<T, 3 ,3 ,T> covTxTyQoP_xyz = .... ;
* ROOT::Math::SVector<C,3>& particleMomentum=...;
* ROOT::Math::SMatrix<T, 4 , 3, T > covPxPyPzE_xyz;
* massOfParticle = ...;
*
* ROOT::Math::SMatrix<T,4,3,T> Jacob =
* JacobdP4dMom (particleMomentum, massOfParticle, chargeOfParticle) ;
* covPxPyPzE_xyz = Jacob * covTxTyQoP_xyz;
*
* @endcode
*
* @param mom (input) the txtyqop vector of track/particle momenta
* @param mass (input) the particle mass
* @param charge (input) the particle charge
* return J (output) the Jacobian for the transformation
* @author Hendrik JAGE hendrik.jage@cern.ch
* @date 2024-04-30
*/
template <class T, class M, class Q>
auto JacobdP4dMom( const ROOT::Math::SVector<T, 3>& mom, const M mass, const Q charge ) {
ROOT::Math::SMatrix<T, 4, 3> J;
auto tx = mom( 0 );
auto ty = mom( 1 );
auto qop = mom( 2 );
using std::abs;
auto p = abs( charge / qop );
auto n2 = 1 + tx * tx + ty * ty;
using std::sqrt;
auto n = sqrt( n2 );
auto n3 = n2 * n;
auto pz = p / n;
auto px = pz * tx; // px = p * tx / n
auto py = pz * ty; // py = p * ty / n
using std::hypot;
auto E = hypot( p, mass );
auto n2_inv = 1 / n2;
auto n3_inv = 1 / n3;
auto qop_inv = 1 / qop;
J( 0, 0 ) = p * ( 1 + ty * ty ) * n3_inv; // dpx/dtx
J( 0, 1 ) = p * tx * -ty * n3_inv; // dpx/dty
J( 0, 2 ) = -px * qop_inv; // dpx/dqop
J( 1, 0 ) = p * ty * -tx * n3_inv; // dpy/dtx
J( 1, 1 ) = p * ( 1 + tx * tx ) * n3_inv; // dpy/dty
J( 1, 2 ) = -py * qop_inv; // dpy/dqop
J( 2, 0 ) = pz * -tx * n2_inv; // dpz/dtx
J( 2, 1 ) = pz * -ty * n2_inv; // dpz/dtx
J( 2, 2 ) = -pz * qop_inv; // dpz/dqop
J( 3, 0 ) = 0.0; // dE/dtx
J( 3, 1 ) = 0.0; // dE/dty
J( 3, 2 ) = p / E * -p * qop_inv; // dE/dqop
return J;
}
// ========================================================================
} // namespace Math
} // end of namespace Gaudi
// ============================================================================
Loading