| 1 |
|
| 2 |
#include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
|
| 3 |
|
| 4 |
|
| 5 |
TwoBodyDecayModel::TwoBodyDecayModel( double mPrimary, double mSecondary ) :
|
| 6 |
thePrimaryMass( mPrimary ), theSecondaryMass( mSecondary ) {}
|
| 7 |
|
| 8 |
|
| 9 |
TwoBodyDecayModel::~TwoBodyDecayModel() {}
|
| 10 |
|
| 11 |
|
| 12 |
AlgebraicMatrix TwoBodyDecayModel::rotationMatrix( double px, double py, double pz )
|
| 13 |
{
|
| 14 |
// compute transverse and absolute momentum
|
| 15 |
double pT2 = px*px + py*py;
|
| 16 |
double p2 = pT2 + pz*pz;
|
| 17 |
double pT = sqrt( pT2 );
|
| 18 |
double p = sqrt( p2 );
|
| 19 |
|
| 20 |
AlgebraicMatrix rotMat( 3, 3 );
|
| 21 |
|
| 22 |
// compute rotation matrix
|
| 23 |
rotMat[0][0] = px*pz/pT/p;
|
| 24 |
rotMat[0][1] = -py/pT;
|
| 25 |
rotMat[0][2] = px/p;
|
| 26 |
|
| 27 |
rotMat[1][0] = py*pz/pT/p;
|
| 28 |
rotMat[1][1] = px/pT;
|
| 29 |
rotMat[1][2] = py/p;
|
| 30 |
|
| 31 |
rotMat[2][0] = -pT/p;
|
| 32 |
rotMat[2][1] = 0.;
|
| 33 |
rotMat[2][2] = pz/p;
|
| 34 |
|
| 35 |
return rotMat;
|
| 36 |
}
|
| 37 |
|
| 38 |
|
| 39 |
AlgebraicMatrix TwoBodyDecayModel::curvilinearToCartesianJacobian( double rho, double theta, double phi, double zMagField )
|
| 40 |
{
|
| 41 |
double q = ( ( rho < 0 ) ? -1. : 1. );
|
| 42 |
double conv = q*zMagField;
|
| 43 |
|
| 44 |
double stheta = sin( theta );
|
| 45 |
double ctheta = cos( theta );
|
| 46 |
double sphi = sin( phi );
|
| 47 |
double cphi = cos( phi );
|
| 48 |
|
| 49 |
AlgebraicMatrix curv2cart( 3, 3 );
|
| 50 |
|
| 51 |
curv2cart[0][0] = -rho*cphi;
|
| 52 |
curv2cart[0][1] = -rho*sphi;
|
| 53 |
curv2cart[0][2] = 0.;
|
| 54 |
|
| 55 |
curv2cart[1][0] = cphi*stheta*ctheta;
|
| 56 |
curv2cart[1][1] = sphi*stheta*ctheta;
|
| 57 |
curv2cart[1][2] = -stheta*stheta;
|
| 58 |
|
| 59 |
curv2cart[2][0] = -sphi;
|
| 60 |
curv2cart[2][1] = cphi;
|
| 61 |
curv2cart[2][2] = 0.;
|
| 62 |
|
| 63 |
curv2cart *= rho/conv;
|
| 64 |
|
| 65 |
return curv2cart;
|
| 66 |
}
|
| 67 |
|
| 68 |
|
| 69 |
AlgebraicMatrix TwoBodyDecayModel::curvilinearToCartesianJacobian( AlgebraicVector curv, double zMagField )
|
| 70 |
{
|
| 71 |
return this->curvilinearToCartesianJacobian( curv[0], curv[1], curv[2], zMagField );
|
| 72 |
}
|
| 73 |
|
| 74 |
|
| 75 |
AlgebraicVector TwoBodyDecayModel::convertCurvilinearToCartesian( AlgebraicVector curv, double zMagField )
|
| 76 |
{
|
| 77 |
double rt = fabs( zMagField/curv[0] );
|
| 78 |
|
| 79 |
AlgebraicVector cart( 3 );
|
| 80 |
cart[0] = rt*cos( curv[2] );
|
| 81 |
cart[1] = rt*sin( curv[2] );
|
| 82 |
cart[2] = rt/tan( curv[1] );
|
| 83 |
|
| 84 |
return cart;
|
| 85 |
}
|
| 86 |
|
| 87 |
|
| 88 |
const std::pair< AlgebraicVector, AlgebraicVector > TwoBodyDecayModel::cartesianSecondaryMomenta( const AlgebraicVector & param )
|
| 89 |
{
|
| 90 |
double px = param[TwoBodyDecayParameters::px];
|
| 91 |
double py = param[TwoBodyDecayParameters::py];
|
| 92 |
double pz = param[TwoBodyDecayParameters::pz];
|
| 93 |
double theta = param[TwoBodyDecayParameters::theta];
|
| 94 |
double phi = param[TwoBodyDecayParameters::phi];
|
| 95 |
|
| 96 |
// compute transverse and absolute momentum
|
| 97 |
double pT2 = px*px + py*py;
|
| 98 |
double p2 = pT2 + pz*pz;
|
| 99 |
double p = sqrt( p2 );
|
| 100 |
|
| 101 |
double sphi = sin( phi );
|
| 102 |
double cphi = cos( phi );
|
| 103 |
double stheta = sin( theta );
|
| 104 |
double ctheta = cos( theta );
|
| 105 |
|
| 106 |
// some constants from kinematics
|
| 107 |
double c1 = 0.5*thePrimaryMass/theSecondaryMass;
|
| 108 |
double c2 = sqrt( c1*c1 - 1. );
|
| 109 |
double c3 = 0.5*c2*ctheta/c1;
|
| 110 |
double c4 = sqrt( p2 + thePrimaryMass*thePrimaryMass );
|
| 111 |
|
| 112 |
// momentum of decay particle 1 in the primary's boosted frame
|
| 113 |
AlgebraicMatrix pplus( 3, 1 );
|
| 114 |
pplus[0][0] = theSecondaryMass*c2*stheta*cphi;
|
| 115 |
pplus[1][0] = theSecondaryMass*c2*stheta*sphi;
|
| 116 |
pplus[2][0] = 0.5*p + c3*c4;
|
| 117 |
|
| 118 |
// momentum of decay particle 2 in the primary's boosted frame
|
| 119 |
AlgebraicMatrix pminus( 3, 1 );
|
| 120 |
pminus[0][0] = -pplus[0][0];
|
| 121 |
pminus[1][0] = -pplus[1][0];
|
| 122 |
pminus[2][0] = 0.5*p - c3*c4;
|
| 123 |
|
| 124 |
AlgebraicMatrix rotMat = rotationMatrix( px, py, pz );
|
| 125 |
|
| 126 |
return std::make_pair( rotMat*pplus, rotMat*pminus );
|
| 127 |
}
|
| 128 |
|
| 129 |
|
| 130 |
const std::pair< AlgebraicVector, AlgebraicVector > TwoBodyDecayModel::cartesianSecondaryMomenta( const TwoBodyDecay & tbd )
|
| 131 |
{
|
| 132 |
return cartesianSecondaryMomenta( tbd.parameters() );
|
| 133 |
}
|
| 134 |
|
| 135 |
const std::pair< AlgebraicVector, AlgebraicVector > TwoBodyDecayModel::cartesianSecondaryMomenta( const TwoBodyDecayParameters & tbdparam )
|
| 136 |
{
|
| 137 |
return cartesianSecondaryMomenta( tbdparam.parameters() );
|
| 138 |
}
|
